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SYNOPSIS 


A 2— D THERMO- MECHANICAL FINITE ELEMENT MODEL FOR RESIDUAL 
STRESS DETERMINATION DURING WELDING AND ANNEALING 
A Thesis Submitted 

In Partial Fulfilment of the Requirements 
for the Degree of 
MASTER OF TECHNOLOGY 

by 

B.K. DUTTA 

Department of Mechanical Engineering 
Indian Institute of Technology, Kanpur 
July, 1983 


A formulation for thermo-plastic analysis of plane 
strain case by finite element technique, has been derived, A 
computer program has been developed for elasto-plastic and 
thermo-plastic analysis of 2.-D structures. 

The formulation for thermo-plastic analysis of plane 
strain case and the computer program have been then modified 
for taking into consideration the effect of creep. 

The program has been tested for number of sample 
problems, including problems involving creep. Program is used 
to. compute stress transients during welding and annealing of 
two long plates. The results have been compared with the 
available results and found in good agreement. 

The program is then used tp compute stress transients 
during welding and annealing of a long cylinder and results 
are reported in the thesis . 
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INTRODUCTION 


1.1 General Review 


1.2 A Problem Definition 


1.3 Computational Task 


1.4 The Finite Element Method 



1 . 1 General Review 


Behaviour of solids under different conditions is 
always a field of interest for engineers and scientists. 
Before the introduction of digital computers/ the application 
of solid mechanics largely rested on the simplif ication of 
analytical formulations/ which are sometimes approximate in 
themselves. The gap between application and theory made the 
whole spectacle unreal. The development of digital computers 
and their application have brought radical changes in this 
area. The literature of last ten years clearly reflects this 
change. This is due to the important role played by numerical 
methods. Today an ever increasing number of engineers and 
scientists are using classical equations of solid mechanics 
on digital computers for solving real-life problems and the 
power of these equations is being extended beyond the dreams 
of their inventors. As a result/ one hears frequently from 
engineers, questions as to the arbitrary idealisation of the 
classical analytical equations . A better appreciation of 
the idealisation involved and ways of inproving them are 
often presented. 

1.2 A Problem Definition 

Modern industries fabricate large components (such 
as calandria, end shield etc. in case of nuclear industry), 
which require different mechanical processes, such as welding, 
casting, grinding etc. These operations in turn introduce 
large residual stresses of different nature. Determination 
of these residual stresses is of utmost interest to designers 



as well as fabricators, as they influence the buckling 
strength of the structure and also structure is susceptible 
to brittle fracture, which gets enhanced in presence of nuclear 
radiations. The fatigue life of the structures is also aff- 
ected considerably and structures are susceptible to cold 
cracking because of high nil! ductility temperature (NDT) if 
large residual stresses are introduced during manufacturing. 

Complexity of the determination of residual stresses, 
are characterized as follows: 

1. Large temperature change in a small area 

2 . Inelastic deformation 

3 . Temperature dependency of material properties 

4 . Phase change of the some portion of the component 

5 . Complex boundary conditions . 

Before the< development of digital computers, most of 
the work for determining residual stresses were confined to 
either experimental or analysis of a very sirtplified one- 
dimensional model which is solvable by conventional mathema- 
tics. But the introduction of digital computers brought a 
revolution in this area and engineers are now in a position 
to compute these stresses more accurately and reliably by 
the help of powerful numerical techniques, such as the 
finite difference and finite element methods. 

1.3 Computational Task 

During mechanical processes, the thermomechanical 
transients control the evolution of the residual stresses and 
deformations. Therefore, the thermal stress analysis of the 



mechanical processes requires both the time hystory analysis 
of temperatures and associated stresses. These are due to: 

i) heat flow to the component and 
ii) quasi- static response behavior of the struct ural component 

For an efficient numerical solution, these two compu- 
tational steps are combined and fully integrated. The two 
field problems are resolved sequentially within the staggered 
solution scheme for transient heat transfer and thermoelastic- 
viscoplastic material behavior on the basis of interpolation 
in space and time. Therefore, the numerical aspects of stab- 
ility and accuracy govern the appropriate selection of time 
steps for the incremental solution of the response history. 

The computational task for determining residual 
stresses reliably is complicated by several facts: 

1. The structural component is subjected to highly 
localised heat flow in the form of conduction, convection and 
radiation. 

2 - The fusion of the molten bead with the base mater- 
ial involves a moving contact problem with phase changes and 
moving boundaries. 

3. The mechanical response during heating and 
cooling behavior spans the entire temperature range from room 
to melting temperature and vice versa. 

4. The change of the internal material structure 
introduces latent heat effects on the thermal side as well as 
additional volume changes and transformation plasticity in 
the region of rapid cooling. 



5. The degradation of stiffness and strength during 
heating causes full unloading and stress relief in the heat 
affected zone of high temperature. 

6. The residual stresses depend on the entire history, 
thus the full heating-cooling cycle has to be traced taking 
due account of the temperature and deformation hystory. 

1.4 The Finite Element Method 

In one operation, it is difficult for a human mind to 
analyse and understand the behavior of its complex surroundings 
and creations. Hence it is a common tendency to divide a 
system to subsystems or 'elements' , for which behavior is 
readily understood. The collective behavior of all the sub- 
systems is an indication of the original system. This natural 
way of analysing a system is adopted by engineers, scientists 
and even economists for many years. 

In general for getting exact behavior of a system, 
one has to go for infinite subdivisions and the problem can 
only be defined using the mathematical fiction of an infini- 
tesimal, which ultimately leads to differential equations or 
equivalent statements. Such systems are called "continuous 
system" . If a finite number of subdivisions can predict the 
exact behavior of the original system, the system is termed 
as "discrete system" . 

With the introduction of high speed digital computer , 
a discrete problem can be solved readily even if the number 
of elements is very large. This gives rise to a tendency of 
approximating a continuous system by a discrete system of 



large number of elements. Since capacity of all computers is 
finite, hence this 'large' cannot be infinite. Hence conti- 
nuous problems can only be solved exactly by mathematical 
manipulation. Here, the available mathematical techniques 
usually limit the possibilities to oversimplified situations . 

The various methods had been proposed time to time 
to approximate a continuous system by a discrete system, also 
termed as 1 discretization* . All these methods involve an 
approximation which, hopefully, is of such a kind that it 
approaches, as closely as desired, the true continuum solution 
as the number of discrete variables increases. 

The discretization of continuum problems has been 
approached differently by mathematician and engineers. The 
first have developed general techniques applicable directly 
to differential equations governing the problem, such as 
finite difference technique. On the other hand, engineers 
often approach the problem more intuitively by creating an 
analogy between real discrete elements and finite portions 
of a continuum domain. 

It is for the engineering 'direct analogy* view that 
the term 'finite element' has been born. Clough appears to 
be first to use this term, which implies in it a direct use 
of standard methodology applicable to discrete systems. Both 
conceptually and from the computational viewpoint, this is 
of the utmost importance. The first allows an improved under- 
standing to be obtained? the second the use of a tonified 
approach to the variety of problems and the development of 
standard computational procedures. 
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CHAPTER 2 

LITERATURE SURVEY AND SCOPE OF THE PRESENT WORK 

2.1 Development of Finite Element Technique 

2.2 Development of Numerical Techniques for Plasticity 
Analysis 

2.2.1 Some Basic Relations 

2.. 2. 2 Solution of Plasticity Problems 

Tangential Stiffness Method 
Initial Strain Method 

Initial Stress Method 

2.3 Development in Thermo-plastic Analysis 

2.4 Development of Residual Stress Analysis During Welding 

2.4.1 One- Dimensional Analysis 

2.4.2 Two-Dimensional Analysis 

Thermal Analysis 

Stress and Distortion Analysis 

2.5 Scope of the Present Work 



2.1 Development of Finite Element Method 


The development of the finite element method (FEM) 
can be described in the following five stages. 

Stage 1 - Initial intuitional development in which Turner ^ ^ ^ 

(1956.) , Argyris^^ (1955)/ Clough*' 3 ^ (1960) presented the 
method as an approximation of continuum in the context of 
structural analysis and Clough for the first time used the 
name " finite element" . 

Stage 2 - Wide recognition jf FEM in the field of structural 

( 4 ) 

mechanics in which the basis for compatible model by Melosh 
(1963)/ equilibrium model by Fraeijs de Verbeke^ 5 "* (1964) and 
hybrid model by Pian^ (1964) depending upon different varia- 
tional principles in elasticity were derived. Numerous finite 
element models and their applications were being reported 

everyday and a proper classification and treatment of the 

( 7 ) 

subject is given in reference by Argyris (1965), Zienkie- 

( 8 ) ( q ) 

wicz and Cheung (1967) and Prezemieniecki 1 (1968). 

Stage 3 - The mathematical interpretation of the FEM and 
convergence are consolidating the fundamentals and are making 
FEM a powerful, universal technique. Ziehkiewicz and his co- 
workers [l] have pioneered many non- structural applications 

of the FEM. Thus other useful processes of finite element 
( 10 11 ) 

derivation 7 exist which widen its scope beyond the 
problems in which a well defined variational functional 
exists. 

Stage 4 - The widening of scope of FEM was first pioneered 

(12) 

in the classic paper by Irons (1966) in which future 



strategy of handling finite element programs is outlined. 

As a result, a family of isoparametric and allied numerically 
integrated elements have been investigated. Most of these 
are reported by Zienkiewicz [ 1 ] . Now many practical aspects 
of reduced and selective integration rules and best sampling 
points for stresses etc. have come to the forefront. The 
research workers have been liberated from the time— wasting 
effort of deducing and programming various matrices afresh 
for each new problem and have directed their attention to 
efficient programs capable of breaking 'cost barriers'. The 
effect of this phase of development on the future applications 
and extensions such as in nonlinear analysis is enormous and 
f ar-r eaching . 

Stage 5 - The developments in FEM today can be seen mainly 

in three directions. In one direction, we find engineers, 
mathematicians and scientists are engaged in improving the 
present solution methods, assembly procedures, stiffness 
computing procedures and also different convergence criteria 
by the help of numerical mathematics. Impact of this work 
can be seen from the development of frontal solution proce- 
dure, skyline assembly technique [37], macro stiffness compu- 
tation procedure, R.L. Taylor [l] and also from the publi- 
cation of enormous number of books in FEM written in purely 
mathematics point of view. 

In the second direction, we find engineers and 
scientists engaged in applying existing finite element solu- 
tion technique to solve engineering and scientific problems. 
Such as use of FEM for analysing temperature distribution in 
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large components , mechanical processes [ll] , large civil 
structures, fatigue [9], creep [18], neutron flux distribution 
in reactor core etc. Even though ASME Cades have not speci- 
fically advocated the use of FEM, it is the most popular tool 
today to compute the detailed stress distribution as required 
by the ASME Code, Section III. 

In the third direction, engineers and scientists are 
engaged in applying FEM in more and more new fields such as 
fluid mechanics, magnets and electro dynamic problems [l] etc. 

2.2 Development of Numerical Techniques for Plasticity 
Analysis 

Theoretical as well as practical interest in plasti- 
city has stimulated much work in the numerical solution of 
elasto-plastic problems by the FEM and other procedures . The 

literature on small strain displacement elastic-plastic solu- 

(13—15) 

tions using FEM is now large [6,33,36], difference lie 

between the apparent forms of the non-linear solution algo- 
rithms used, and the form of constitutive relations postulated. 

Before 1971, in most of the literature on the subject, 
simple, constant stress, elements have been used. This was 
largely motivated by the situation in which an exact integ- 
ration of an element containing both elastic and plastic 
regimes would be difficult or impossible. Such difficulties 
disappear in the case of numerically integrated elements. 
Numerical integration was then introduced extensively into 
the formulation of complex elements [l] in linear elastic 
analysis and the higher order elements were proved to be 
remarkably efficient. Adoption of isoparametric formulation 



for materially nonlinear problems was done for the first time 
by G .C . Nayak (1971) [6] and showed that even greater advan- 
tages arise once a smooth distribution of the displacement 
throughout the element is postulated by a shape function 
definition. So state of stress at each integrating point can 
be examined individually and the resulting element properties 
obtained by a suitable weighted summation over the integrating 
points . 

2.2.1 Some Basic Relations 

To discuss developments in solution procedures for 
plasticity analysis , some basic relations of plasticity, are 
given here for further reference. Derivations of these 
relations are given in Chapter 3 . 

incremental strain-displacement relation £de } = [b] {d5} 







(2 

.1) 

incremental equation of equilibrium 

CdR} 

- / [b] T 

(dq 

dV 





V 









(2 

.2) 

tangent stiffness matrix [Krp] = 

/ 

V 

[Bf 

Mep !>] 

dV 

(2 

.3) 

incremental plastic strain {de}^ 

= 

deP 

, 3 0 

* 3 o' ^ = 

d£ 

H' 

? il. 

30 

} 



9 a 

i 30 i 






* go 


(2 

.4) 


= 


H* 




incremental elastic-plastic stress strain relation 
{doj = [D] ep Sde; = [ [D] - [D] p ] Jde ; 

= (do } e - {dd} p (2.5) 


where 
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l>J 5-tl-i L#J [D] 


FI 2 + 


a o 5 u 3 cr- 

iH-J m 


{iUL i 

30 3 


(2.6) 


Hence any plasticity problem now involves the solution 
of equilibrium Equation (2.2) together with incremental 
stress strain relations given by Equations (2.5) and (2.6). 

2.2.2 Solu tion of Plastic ity J? roblerns 

Tange ntia l Stiffness Met hod [33] - By substituting { daj and 

{ds} from Equation (2.5) and Equation (2.1) into equilibrium 
Equations (2.2) , we have 


idR} - / [B] T [D] [B] dV {d 5 } = {0} 

v 

which from Equation (2.3) reduces to 


(2.7) 


{dR} - [K.J.] Jdo} = co; 


( 2 . 8 ) 


hence 


id6} 


[kJ - 1 idRi 


(2.9) 


This method of solving was presented by Pope v } 

(1965) for plane stress problems for strain hardening von 
Mises solid. He advocated the use of small increments of load 
which just caused yield in the next simple triangular element . 
Yamada, Yoshimura and Sakurai (1968) [32] used the method 
presented by Pope by allowing one triangular element to yield 
at a time. They have clearly shown that the problem is very 
much dependent on the mesh type and sizes and there can be 
some elastic elements surrounded by plastic ones near to the 


elastic-plastic interface 
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Marcal and King (1967) [33] advocated a large size 
increment of load (based on initial yield load factor L ) and 
solved the Equations (2.9) by using a mean value of tangen- 
tial stiffness (partial stiffness) for elastic perfectly 
plastic and strain hardening Von Mises solid with simple 
triangular elements in plane stress , plane strain and axy— 
symmetry. The necessity of partial stiffness coefficients 
arose because of the 'transition region' which is usually 
adjacent to the elastic plastic interface and for correct 
solution of Equations (2,9). The transition region must be 
correctly known. 

Partial stiffness matrix for an element under tran- 
sition is calculated by 

C D ep ] = O] - (1 - m) [D] p (2.10) 

0-0 

where m = r p - - • and a > a > a „ (2.11) 

a n “ a n-l n y n " 1 v * J 

More explanations on this method can be seen in Chapter 3 . 


in i ti al Strain Method - Another modification of Equations 
(2.2) is based on the idea of modifying the right hand side 
of elastic equations of equilibrium to compensate for the 
fact that the plastic strains do not cause any change in 
stress and can be treated in the same way as the thermal 
strains . Singe 

{A a } = [D]{Ae} e = [d] ( ( A £ } fc -CA£} P ) (2.12) 


hence using Equation (2.4) in Equation (2.12), we have 


{A<?} = [D] ( {A e } 


iM i {AO} 
,t 1 UoJ i 


H' 


) 


(2.13) 
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Using Equation (2.13) in Equation (2.2) we have 

{ARJ = [ K 0 ] { A<3 } ~ / [Bf [D] (~- L-—J {AO}) dV 

V o 

(2.14) 

where 

[*y = i [b] t [d] [bj dv 
v 

Hence 


[KJ {Ao} = {AR} + {AP } (2.15) 

where 

JAP} = / ~ ([b] T [d] } L~-J {AO}) dv (2.16) 

Obviously solution of Equations (2.15) is impossible 
if H' =0 (that is, no hardening). 

Gallagher et al^^ (1962) and Argyris^®^ (1965) have 
effectively used Equation (2.16). The essential difference 
from the tangential stiffness method is the evaluation of {AP } 
by iterations. 

Evidently the method satisfies the equilibrium equa- 
tions and stress-strain relationships. The method is not 
applicable for ideal elastic-plastic material but can be 
made to work with hardening non-associated flow laws. For 
solution, [ K ] ^ can be stored once and for all. It is, 
therefore, not likely to be as costly as the tangential 
stiffness method. However as the yielding of material 
increases the number of iterations required are more. 

Initial stress Method - Since the 1 initial strain method' 
is not applicable to ideally elastic-plastic material an 
alternative method called, 1 initial stress method' was 



(1968- 


introduced by ZIenkiewicz, Vaf liappan and King ^ 19 ^ 

1969} . Here Equations (2.2) are written as 

rii ; =* [K] iAo} - f [b] T {AO} p av (2.17) 

u V 

where 


{A (■ -- [d] p {A©} 

hence 'he .initial load 


(2.18) 


{A?} = [l<] P {Ao } --- f [ B] T t AG} P dV (2.19) 

V 

This method has the advantage of being easily inter- 
preted physically and is simple and presents no problems for 
cyclic loading as the unloading proceeds on purely elastic 
basis. There are no problems as regards the consideration 
of nor-associated flow laws and has an obvious advantage of 
operating with [k q ] ^ stored once and for all. However its 
preference over other methods such as tangential stiffness 
method depends mainly on the rate of convergence. The initial 
stress method converges very rapidly in the beginning and 
becomes very slow near general yield. 

There are various other methods which are being used 
such as method of residual forces [ 6] , alpha constant stiff- 
ness method " , method of plastic multiplier [6] etc. But 

the three methods explained above seems to be popular . A 
computer program developed by the present author, uses the 
tangential stiffness method as devised by Marc'll and King. 

One of the reasons of using this method is variations in 
material properties with loading (because of their temper- 
ature dependency) and hence matrix [K q ] < ” 1 cannot be stored 



and used for" further iterations- Since it is now necessary 
to calculate stiffness matrix for every loading step, it is 
better to calculate [k] t instead of [k] , as in former case 
convergence is faster. 

2 . 3 Development in Thermo-Plastic Analysis 

A thermo-plastic analysis can be treated very much 
similar to elasto-plastic analysis . But in case of thermo- 
plastic analysis beside elastic and plastic strains, two more 
strains corresponding to the temperature change and the 
temperature dependent material properties, are to be consi- 
dered. Tangential stiffness method seems to be a popular 
choice because of necessity of computation of stiffness matrix 
in every loading step due to the dependency of the material 
properties on temperature. One of the oldest paper in this 
analysis is that of UEDA et al (1971) [19]. The formulation 
presented in [19 j and also in [l7] by the same authors, seems 
to be applicable for 1-D, 2-D plane stress and axisymmetric 
and 3-D structures. But no formulation was reported for 
plane strain case in which a systematic elimination of z 
direction incremental stress is required. Subsequently many 
authors used the above formulations in their work, such as [ 9] , 
[12] , [3] etc. For plane strain case a method is suggested 
by Friedman (1975) [20], by postulating an expression for 
d&P . In the present work systematic derivation of different 
expressions for thermo-plastic analysis of plane strain case 
is done and shown in Chapter 3 . 



2*4 Development, of Residual Stress Analysis During Welding 

The welding process is one of the roost widely used 
structural fabrication techniques, and has therefore received 
much study as a mechanical process. The objective of such 
work, is the identification of problems associated with welding . 
Among these the residual stresses and distortions resulting 
from the extreme thermal conditions during welding are one 
most important. The residual state of stress has a strong 
influence on the strength and life of the resulting structure, 
and welds are usually the critical determinants of the overall 
structural integrity. Both experimentally and analytically, 
the problem poses severe difficulties, and inspite of the 
extent of the experimental effort, there have been few analy- 
tical studies of the residual stresses in a weld. Analyti- 
cally most of the studies available can be divided into two 
categories. These are the studies done before the development 
of FEM and the studies done after the development of FEM. 

Before the development of FEM most of the available work is 
one-dimensional and uses a line solution for the thermal 
analysis. Whereas after 1970s, work done in this context is 
of 2-D analysis by FEM. Studies on the transient thermal 

stresses during welding started in the 1930s. Boulton and 
( 21 ) 

Lance-Martin discussed the transient thermal stresses 
that occur along the edge of a plate during welding. 

The first significant attempt to use a computer in 
the analysis of thermal stresses during welding was done by 
Tall ( 22 ' 23 ) j_ n a ph.D. thesis in 1961. He developed a simple 
program on thermal stresses during bead welding along the 
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center line of a strip. The temperature distribution was 
treated as two-dimensional; however, in analyzing stresses 
it was assumed that (1) longitudinal stress, d is a function 
of the lateral distance Y only and (2) that c and t are 

x jDl 

zero. 

( 24 ) 

In 1968 Masgbuchi et al v ' of Battelle Memorial 
Institute developed based upon Tail's analysis, a FORTRAN 
program on the one-dimensional analysis of thermal stresses 
during welding. Since 1970 the computer analysis of transi- 
ent thermal stresses during welding has become more common at 
several laboratories around the world [ 4 ]. 

Commission X (Residual stress, stress relieving and 
brittle fracture) of the International Institute of Welding 
established in 1972 a working group of "Numerical analyses of 
stresses, strains and other effects produced by welding" . 

The working group has prepared reports covering the studies 

( 25 ) 

made in various laboratories of the World . 

In this section a general one- dimensional analysis 
done before 1970s and two*> dimensional analysis after 1970s 
and inherent assumptions in both the analysis separately are 
presented. 



Temperature distribution around the moving arc is 
calculated as a two-dimensional heat-conduction problem. The 
temperature distribution can be expressed by the following 
equation [ 4] ; 



2n\ 
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where 

\H; = temperature change 
q = intensity of the linear heat source 
X = heat conductivity 
X = thermal diffusivity 
V = arc travel speed 

= x - Vt - longitudinal ordinate in the moving 
coordinate system 
= time 

;? = V g 2 + Y 2 

K 0 (2) - zero order modified Bessel function of the 

second kind „ 

The Intensity of the linear heat source, q, is expre- 
ssed by: 

r) V I 

q = 

where r? ~ thermal efficiency of welding arc, or arc 
efficiency 
V = arc voltage 
I = arc current 
T ~ plate thickness 

For calculating stresses, the field is divided into 
a set of transverse strips of width h Q , as shown in Figure 
2.1* The time intervals represented by the strip width must 
be narrow enough so that the temperature and thermal stress 
for each increment may be regarded as being constant. Since 
the greatest changes In temperature occur near the arc, 
narrow strips are used in areas near the arc. 
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Key assumptions in the one- dimensional programs are 

that 

(a) both sides of each strip remain straight during welding, 
which produces longitudinal stresses, o . creating a sort 
of plane- strain condition in each strip; and 

(b) the longitudinal stresses along the entire length of 
each strip are balanced. 

Figure 2 .2 shows schematically changes of temperature 
and stresses during welding as predicted by 1-D analysis . A 
bead-on-plate weld is being made along the X-axis. The welding 
arc, which is moving at a speed v, is presently located at 
the origin 'O', as shown in Figure 2.2(a). Figure 2.2(b) shows 
temperature distribution along several cross sections. Along 
section A-A, which is ahead of the welding arc, the temper- 
ature change due to welding, AT, is almost zero. Along 
section B— B which crosses the welding arc, the temperature 
distribution is very steep. Along section C-C which is some 
distance behind the welding arc, the distribution of temper- 
ature change is as shown in Figure 2.2(b-2). Along section 
D-D, which is very far from the welding arc, the temperature 
change due to welding again diminishes. Figure 2.2(c) shows 
the distribution of stresses along these sections in the 
x direction. 

Along section A-A, thermal stresses due to welding 
are almost zero. The stress distribution along section B-B 
is shown in Figure 2.2(c-l). Stresses in regions underneath 
the welding arc are close to zero, because molten metal does 
not support loads. Stresses in regions some distance away 



from the arc are compressive, because the expansion of these 
areas is restrained by surrounding metal that is at lower 
temperatures . Since the temperatures of these areas are quite 
high and the yield strength of the material is low, stresses 
in these areas are as high as the yield strength of material 
at corresponding temperatures. The magnitude of the compre- 
ssive stress passes through a maximum with increasing distance 
from the weld or with decreasing temperature. However, 
stresses in areas away from the weld are tensile, and balance 
with compressive stresses in areas near the weld. In other 
words / dy =0 across section B-B , Thus, the stress 
distribution along section B-B is as shown in Figure 2.2{c-l). 
Stresses are distributed along section C-C as shown in 
Figure 2.2(c-2)« Since the weld-metal and base-metal regions 
near the weld have cooled, they begin to shrink, causing 
tensile stresses in regions close to the weld. As the dis- 
tance from the weld increases, the stresses first change to 
compressive and then become tensile. Figure 2.2(c-3) shows 
the stress distribution along section D-D. High tensile 
stresses are produced in regions near the weld, while 
compressive stresses are produced in regions away from the 
weld. The, distribution of residual stresses that remain 
after the welding is completed, is shown in figure. 

2.4.2 Two-Dimensional Analysis 

A two-dimensional model was first developed and 

/ O tL ) 

applied by Hibbitt and Marcel successfully . Comparison 
with available experimental work was made, with gross 
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quantitative agreement . The same model wqs used by Nickel 1 

and Hibbitt [39]/ who compared temperature results (including 

the prediction of fusion zone) with a well instr lamented 

experiment on a complex geometry with good accuracy. A 

similar model has been developed by Friedman [20], who studied 

butt welds, but no experimental comparison has been made. 

(27 ) 

Buyukozturk and Hibbit applied this technique to a large, 

multipass butt weld, again without experimental verification. 
In all these cases the authors reported that their main diff- 
iculty and source of possible error was the lack of basic 
material and parametric data. 

Thermal Analysis 

Heat Input - The mechanisms involved in heat input 
are extremely conplex and physics of the weld arc are not 
fully understood [ll] . The heat of the arc may be related 
to 

q = n v i 

where Q is the net heat flux, Tj is the arc efficiency, V is 
the electrical input voltage and I is the electrical input 
current . 

The distribution of flux on the facing surface is 
assumed to be normal and radially symmetric. That is 

2 

— c r 

q = q o e 

where q Q and c are chosen based on the radius of the welding 
electrode. The most effective method of obtaining these two 
parameters is from a well instrumented test. 



In most welds it is readily shown that the arc speed 
is large compared to heat transfer rates, so that the heat 
transfer parallel to the line of motion of the electrode may 
be neglected. That is, the temperature history of any section 
normal to the line of motion of the weld head is the same as 
that of any other such section, except for a shift in time. 
This is an important consideration, since it allows a reduc- 
tion from a three-dimensional to a two-dimensional thermal 
analysis . Hence 


2 -,2 , ^ .2 

t . \ -CS -CV (t — to) 

q(s, t) = q o e e ° 

where ‘s' measures distance, on the surface, from the center- 
line of the weld electrode. V is the electrode velocity and 
t is the time of peak heating for this section. 

Heat Transfer - The mechanism which carries the 
heat applied to the surfaces described above, into the bulk 
material of the filler and base material adjacent to the 
joint is that of conductive heat transfer. The Fourier law 
q 3 = -K . .T . . governs this mechanism, where K, .(T) is the 
conductivity . 

In addition, specific and latent heat effects must 
be included, since the problem is a transient one, with 
melting and re-solidification occurring in the material. 

Since both the specific heat and conductivity are known to 
vary with terrperature, it is necessary to have their values 
available for the complete temperature range of the analysis. 
Usually the investigators resort to extrapolation of known 



data,, since high temperature values are difficult to measure 
experimentally [19,20]. 

The FEM as conventionally applied to heat transfer 
problems is not well suited to latent heat calculations. The 
reason for this is that the method is based on the assumption 
of a local solution, T = (x) T N , which is assumed C 1 conti- 

nuous on the interior of the elements . When phase change 
take; place at a unique temperature (as in a pure phase) the 
latent heat requires a discontinuity in the temperature grad- 
ients » Such discontinuities are only admitted at the element 
interfaces, so that it is not realistic to resolve the fusion 
front any finer than the element dimension. Apparent incr- 
ease in heat capacity seems to be a popular method for phase 
change consideration [l4] . 

Boundary Conditions - The third mechanism of impor- 
tance in the heat transfer is the surface boundary conditions. 
Two heat transfer mechanisms are usually considered: radia- 
tion effects dominate at higher tenperatures , but later in 
the cool-down history the primary mechanism is convection. 

The difficulty in the latter case is the inside variation in 
convection coefficient depending on surface condition and 
air flow. 

Galepkin finite element method [l] and central 
difference (Crank-Nicholsen) time stepping algorithm [l3] 
seem to be again popular choice [11] . 

St ress and Distortion Analysis - The stress analysis also 
utilizes a finite element geometric model and piecewise 



linear time-stepping. The input for this phase of the anal- 
ysis is the temperature history defined throughout the joint 
region by the previous thermal analysis . 

Finite Element Model - In most cases the geometric 
model must be modified for the stress analysis [20], This is 
because the thermal analysis is a very local solution. But 
in order to predict residual stress and distortion accurately, 
it is important to model the flexibility of the rest of the 
structure, since the constraint provided by the structure 
outside the region in the immediate neighborhood of the joint 
plays a significant part in the growth of residual stress. 

Material Properties - The residual stresses res ul t 
from inelastic response of the material to the extremely 
severe thermal history of the joint constitutive character- 
ization of the material is difficult since all temperature 
from room temperature to well above the melting point are 
present. In all the work to date, the assumption is made 
that the material may be treated as a time- independent elastic- 
plastic material during welding, with temperature dependent 

yield and elastic moduli. This characterization was justi- 
( 26 ) 

fied in v by the argument that welding times are short so 
that any time- dependent behavior is restricted to sharp 
primary creep effects of extremely short characteristic 
times, which cannot be distinguished easily from yield in 

m 

a material characterization experiment. 

An interesting aspect of the weld process is the 
presence of sharp temporal temperature gradients. The mater- 
ial cools so rapidly that it is clear that unstable material 



phases must be present, so that the material properties in 

regions, where this happens, must be quite different from 

those of the bulk material. This point has never received 

attention in the welding analysis that have been performed, 
(27 ) 

except in , where the authors use a different material 
characterization in the heat-affected zone, based on post- 
weld surface hardness measurements. 

Elastic-Plastic Solution Procedure - The use of 

the FEM for non-isothermal elastic-plastic analysis is by 

now a routine procedure, and the weld problem is merely a 

rather severe example of such an analysis . All the authors 

methods 

cited use tangent modulus^[7 ,11 ,12,19,20] where by the system 
stiffness matrix is reformed at each step in order to account 
for active plastic loading, or elastic unloading from such a 
state. The formulation uses a piecewise linear approximation 
for the classical equations of plasticity, and some scheme 
is always used to process the thermal analysis results to 
ensure sufficiently small steps in the stress analysis to 
satisfy this linearization requirement. 

2.5 Scope of the Present Work 

Welding of long plates, cylinders etc. can be 
treated as a plane strain case for stress transients anal- 
ysis. In the present work, a complete derivations of diff- 
erent expressions for thermo -plastic analysis of plane strain 
case by FEM, has been carried out . Incremental stress in 
Z direction is eliminated during this derivation, because of 
its dependency on other inplane stresses. A computer program 



THESIS has been developed for thermo-plastic analysis of 
plane- stress , plane- strain and axi symmetric structures depen- 
ding upon the formulations derived. The same program has 
been also modified for taking into consideration the melting 
of part of the structure. The program has been tested for 
different sample problems and then used to compute stress 
transients during butt welding of two long plates and longi- 
tudinal welding of a cylinder. 

Upper derivation of different expressions for thermo- 
plastic analysis of plane strain case has been repeated for 
taking into consideration creep. Pour creep laws such as 
power creep, experimental creep, time hardening creep and 
strain hardening creep laws, have been used in the formula- 
tion for better representation of creep phenomena of any 
materials during any temperature range. Program THESIS has 
been then modified for thermo-plastic analysis with creep 
of 2-D structures. The program has been again tested for 
different sample problems and then used for computing stress 
transients during annealing of two butt welded plates and 
longitudinal welded cylinder. 

Eight noded isoparametric element has been used for 
the finite element modelling of the structure. 2x2 Gauss 
point integration, Gaussian elimination solution procedure 
and macro stiffness computing procedure for fast computation 
of element stiffness matrices, as suggested in [l] , have been 
used in the program THESIS. Tangential stiffness solution 
procedure for plasticity analysis has been used as discussed 
in Section (2.2.2). 



CHAPTER 3 


ELASTO-PLASTIC AND THERMO-PLASTIC ANALYSIS OF 2-D STRUCTURES 
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3.2.1a Tresca-Yield Criteria 

3.2.1b Von-Mises Criteria 
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Yield Criteria 
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Thermo-Plastic Analysis of Plane Stress and 
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3.7/ Thermo-Plastic Analysis of Plane Strain Problems 

3.8 An Iterative Solution Algorithm 

3.8.1 A Solution Step 

3*8.2 Elastic- Plastic Transition 

3.8.3 Calculation of Unbalance Force 

3.9 A Computer Program and Check-Out Problems 

3*9.1 Thick Cylinder Under Internal Pressure 
(A Plane Strain Case) 

3.9.2 Thick Cylinder Under Radial Temperature 
Gradient (Plane Stress as well as 
Axisymmetric Case) 

3.9.3 Thick Cylinder Under Radial Temperature 
Gradient (A Plane Strain Case). 



3.1 Introduction 


In the present chapter, four sets of formulations are 
presented for elasto-plastic and thermo-plastic analysis of 
2-D structures. The main aim of these formulations is to 
find egressions for elasto-plastic material matrix [d ] 
and corresponding load vector. These expressions are used 
in finite element formulation presented in Appendix 1, for 
ultimate analysis of entire structure. In Section 3.2 basic 
relations for plasticity are presented for further reference. 
In Section 3.3, a general approach to finite element for 
elasto-plastic and thermo-plastic analysis is discussed. A 
formulation is presented in Section 3.4 for elasto-plastic 
analysis of plane stress and axisymmetric structures. In 
Section 3.5, above formulation is modified for plane strain 
case. Formulation for thermo-plastic analysis of plane stress 
and axisymmetric structures is presented in Section 3.6, which 
is an extension of Section 3.4. Finally a new approach for 
thermo-plastic analysis of plane strain case is presented 
in Section 3.7, developed during the course of this work. 

An incremental solution algorithm is given in Section 
3.8. Section 3.9 deals with solved examples to show capa- 
bility of these formulations to attack different types of 
problems . 

3.2 General Theory of Plasticity 

•Plastic' behavior of solids is characterized by a 
non-unique stress-strain relationship, as opposed to that 
of linear elasticity. One definition of plasticity may be 
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the presence of irrecoverable strains on load removal. If 
uniaxial behavior of a material is considered, as shown in 
Figure 3.1(a), a non-linear relationship on loading alone 
does not determine whether non-linear elastic or plastic 
behavior is exhibited. Unloading will immediately discover 
the difference with the elastic material following the same 
path and the plastic material showing a history dependent/ 
different path. Many materials show an ideal plastic behavior 
in which a limiting yield stress, d , exists at which the 
strains are indeterminate, as shown in Figure 3.1(b). For 

I 

all stresses below a linear (or non-linear) elasticity 

relationship is assumed. A further refinement of this model 

is one of hardening/softening plastic material in which the 

yield stress depends on some parameter, such as s as shown 

P 

in Figure 3.1(c). 

Because of non-uniqueness of stress- strain relation- 
ship, one has to consider the following three points, before 
proceeding to plasticity analysis . 

3.2.1 Yield Criteria 

A yield criteria is a hypothesis concerning the limit 
of elasticity under any possible combination of stresses. 

The suitability of any proposed yield criteria must be 
checked by experiment . 

If a point in a ductile material , is subjected to 

the principal stresses , d^, o and is represented by a 

Mohr circle, if the principal stresses are changed to 

( CT . + a _) ( a 0 + d ) and (o- + a ) , then the size of Mohr 
1 m. 2 m 3 m, 



circle, for this new state of stress, remains same but is 


shifted by a distance along the a- axis. The additional 
stresses make up a hydrostatic (tensile or corrpressive) 
stress system. It is found that the absolute size of the 
Mohr circle alone determines the limit of elasticity and is 
independent of its position. This is to say that the yield 
criteria is a function of (c^- o^) , (cr - o^) and (a - a ) 
and is independent of hydrostatic stress component ( 

+ & 3 )/3 . 

Thus yielding occurs when some scalar function of 
the principal stress differences reaches a critical magnitude, 
or to say mathematically: 

f(0 l" °2* °2 “ °3' °3 “ °1 ) = Constant (3.1) 

3.2.1a Tresca- Yield Criteria 

Perhaps the simplest function imaginable which satis- 
fies Equation (3.1) is of the form 

O'. - a. = Constant (3.2) 

i J 

Equation (3.2) interprets when largest of the three 
magnitudes ° 1 ~ ° 2 ' °2 ~ °3 ' °3 ~ °1 attains a critical 

constant value (for a given material) yielding begins. 

Equation (3.2) was first suggested by Tresca in 1864. 
Hence this states: 'yielding occurs when the greatest abso- 
lute value of any one of the three maximum shear stresses 
in the material reaches a certain value' . 

Since yield criteria should be valid for any 
combination of stresses, the constant in Equation (3.2) may 
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be determined from uniaxial stress state. Denoting yield 

shear stress by k for a state of pure shear and tensile yield 

stress by a Yp , since for pure shear stress - a = k 

and the intermediate principal stress cr = 0, the value of 

the constant is 2k. Also for pure tension we have o = u 

1 YP 

and ° 2 = ° 3 = 0, hence for ^ > one gets for yield _ 

ing by Tresca criteria 

!0 1 “ °3 I = a Y P = 2k (3.3) 

3.2.1b Von-Mises Criteria 

Another admissible function, which satisfies Equa- 
tion (3,1) is of the form 

^°l“ °2^ + ( °2 ~ + ( °3 “ ° = Constant (3.4) 

The incorporation of (t^- a ) etc., as quadratic 
terms, dispenses with the need to use the modulus sign. In 
this type of function, each of the principal stresses cont- 
ributes to yielding. This function was proposed by Huber 
(1904), Von-Mises (1913) and by J.c. Maxwell in a letter to 
•^- e l v -in in 1856 £31^j . The criteria was then interpreted by 
Hencky as ‘yielding begins when the shear strain energy 
reached a critical value' . 

The value of the constant in Equation (3.4) can be 
again determined from the simple tension yielding, i.e. 

°l ~ °yp' °2 = ~ 0 and pure shear yielding i.e. a = 

~ ~ = 0. Hence yield criteria becomes 

(0 l _O 2 >2 + ( ° 2 - C 3 )2 + (° 3 - q) 2 - 2°r P = » 2 (3.5) 
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Hence the relation between 0 Yp and k, as per Tresca 
criterion is k = 0.5 <? Yp/ whereas as per Von-Mises criteria 

wF> = 0.577350. 

3 YP 

For a general state of stress ( o, ■ . a . a . r , 

* y z' xy 

T y 2 ‘' T ' Z y^ * tiie Von-Mises criteria becomes 


V2 


[ c 0. 


”/ + <°y' ° 2 ) 2 + <° s - V 2 


+ 6 (T 2 +1 2 + T 2 n 1 / 2 = 

xy yz zx J 


YP 


(3.6) 


The formulations subsequently presented in this 
chapter and in the following- chapters, make use of this 
criteria for yielding as it is the best available yield 
criteria and this includes all the principal stresses in 
one expression. 


3.2.1c Stress Space Representations of Yield Criteria 

In Figure 3.2(a) three mutually perpendicular axes, 

0 0 and 0 g^ are shown. If the principal stresses 

at a point in a body are (c , a , a ), this state of stress 
system is represented by a point P in this ‘ stress space' , 
the coordinates of P being g^, a . The state may be 

written as the sum of the three vectors 0P^ (= o ^) / p^m (= g^) 
and MP (= g^) . 

Let OH be a line, which is equally inclined to all 
the three axes (i.e. direction cosines of OH are 1/V3, 1/V3, 

1/V3). Projection of OP on OH is ON. Hence ON = g^/V3 + 
g 2 /V3+g 3 /V3 and PN 2 = | [ ( ^ _ o 2 ) 2 + (G 2 - a 3 ) 2 

+ ( 0 3 - 



Since according to Von— Mises criteria, yielding 

occurs when (c - a ) 2 + (a - a ) 2 + (o _ a ) 2 = 2 C 2 
x z 2 3 3 1 YP 

(from Equation (3.6)) hence yield locus may be represented 
by a circle of radius PN = sHJ3 . It follows that the 
Von-Mises yield criteria may be represented in principal 
stress space by a circular cylinder of radius C Yp V2/T, 
whose axis is the line through the origin equally inclined 
to the axes of coordinates, as shown in Figure 3.2(b). 
Yielding is unaffected by the magnitude of ON which repre- 
sents a hydrostatic state of stress. ON is the vector sum 
of the spherical components and NP is the vector sum of the 
deviatoric components of the stress. 

The Tresca criteria may be written io - o i = o 

1 3 1 1 YP 

(from Equation (3.3)) and when drawn in principal stress 
space it is represented by a regular hexagonal cylinder which 
is conventionally inscribed within the Von-Mises cylinder, 
as snown in Figure 3.2(b). This can be also demonstrated by 
considering the projection of the stress states on to a 
plane through the origin and perpendicular to OH. This plane 
is called the jr plane and in the original principal stress 
space has the equation + c 2 + = 0. The positive 

principal axes 0 g-^, 0 0 appear in the n -plane, as 

shqwn in Figure 3.2(c), inclined to 120° to each other. 

When the negative axes of principal stress are included, 
the whole area or plane is divided into six equal sectors. 

The representation of Von-Mises criteria and Tresca yield 
criteria on it plane, are shown in Figure 3.2(c). 



3 « 2 » 2 Stress— Strain Relations (Flow Rules) 
3.2 «2a Levy-Rises Equations 


Saint Venant (1870), first proposed that the princi- 
pal ex.es of strain increment coincided with the axes of prin- 
cipal stress [31]. The general relationship between strain 
mcr 3 :nt and the deviatoiric stresses was first introduced 
by L -- (1871) and independently by Von Mises (1913). These 
equa ; ons are now known as the Levy-Mises equations and are 


written as 





a e_ 


zx 


T 

zx 



a\ 


(3.7) 


Since in Equation (3.7/) elastic strains are not 


taken into account, the Levy-Mises relations obviously cannot 
be used to obtain information about ' elastic spring-back* or 
residual stresses, but are useful for unrestricted plastic 
flow. 


3.2.2b Prandtl-Reuss Equations 

The stress-strain relations for an elastic-perf ectly 
plastic solid were first proposed by Prandtl (1924) for the 
case of plane- strain deformation. The general form of the 
equations was given by Reuss (1930) [31]. Prandtl* s equations 

are an extension of the Levy-Mises equations. 

Reuss assumed that the plastic strain increment at 
any instant, is proportional to the instantaneous stress 
deviation and the shear stresses, thus 



41 


d eP 

_ x 

s 

X 


or 


d eP 
— 3C 
s 

y 


d eP 


d sf? • d e p 
... yz _ zx 


d s’ 


_*y 


T 


y* 


'zx 


T xy 


= dX 


d sP , = s . , dx 

ij ij A 


( 3 . 8 ) 


dX i= an i^otantaneous non negative constant of proportion- 
ality hie, may vary throughout loading cycle. 

iho equat ions state that a small increment of plastic 
straj depends on the current deviatoric stress, not on the 
stress increment which is required to bring it about. Also, 
the principal axes of stress and plastic strain increment 
coincide. The equation is only a statement about the ratio 
of the plastic strain increments in the various X, Y, z- 
direations; it gives no direct information about their abso- 
lute magnitude. 


3 . 2.3 


5 ~ ! - n Hardening Theories 


In case of strain hardening plastic materials, though 
elastn^ e iavicr and hence its initial yield condition are 
same as that of a perfectly plastic material, but once plastic 
flow begins, it.s continuation depends on the degree to which 
the material strain-hardens . The criteria for this contin- 
uation may then be formally stated as 


f c (J 2 , J 3 ) 


Constant 


q > o 


where the constant now depends upon the condition of current 
yielding and, therefore, the plastic- stress history. Here 
and are the second and third invariants of the devia- 
toric stress . 





^or the development of any theory of strain hard- 
ening, a basic requirement is an errpirical relation between 
the degree of strain hardening (represented by the current 
yield condition) and the plastic-strain rate. Available 
experimental evidence merely points up the difficulties 
encountered in the formulation of such a relation. The lack 
of significant test data has given rise to some theoretical 
formulations. These theories are as follows. 

3.2.3a Isotropic Strain Hardening Theory 

This theory is due to Hill and is based on the hypo- 
thesis that the degree of strain hardening is a function only 
of the total plastic work and is independent of the strain 
path taken to reach a given stress. Hence this theory 
assumes that the yield surface expands uniformly about the 
origin while maintaining its shape and orientation, as shown 
graphically in Figure 3.3(a). In other words, even if 
yielding occurs in tension, the effect of strain hardening 
on the yield stress in compression is equal to that in 
tension. Thus, tensile strain hardening also numerically 
increases the compressive yield stress and the elastic range. 

3.2.3b Kinematic Strain Hardening Theory 

Another theory of strain hardening by Prager is 
based on the use of a kinematic model. The basic assumption 
of this theory is that, when the material strain-hardens in 
tension, the elastic unloading range is double the initial 
tensile yield stress. Thus, tensile strain hardening 





A conse- 


numerically decreases the compressive yield stress, 
quence of this theory is that the yield surface is free to 
translate as a rigid body while, maintaining its initial 
shape and size, as shown in Figure 3.3(b). 

Yet another theory based on the physical concept of 
grain slip has been developed by Batdorf and Budiansky. 

This theory leads to flow laws which are mathematically diff- 
icult to handle and a few problem solutions that exist are 
not in particularly good agreement with experimental results 
[34] . Other theories combining isotropic hardening and 
kinematic hardening have also been proposed . 

A computer program THESIS, developed during the 
course of this work, described latter in this chapter, makes 
use of isotropic strain hardening theory, i.e. same shape of 
stress— strain curve whether the uniaxial stress is tensile 
or compressive. 

Figure 3.4(a) shows a typical one dimensional stress- 
strain curve for isotropic strain hardening material, the slope 
at any point in this curve is given by 


In Figure 3.4(b), same curve is plotted between one 
dimensional effective plastic strain. Slope at 

any point on this curve is given by 


H* 



(3.9) 


Also from geometry, de^ 


de 


11 

E 



Using the above three relations, we have 



(3.10) 


3 * 3 ~-j^- n:Lte El ement Approach to Thermo- Mechanical Problems 

In this chapter from Section 3.4 onwards dif ferent 
formulations are presented for different cases of 2-D struc- 
tures for elasto-plastic and thermo-plastic analysis. Von- 
Mises yield criteria, Prandtl-Reuss flow rule and isotropic 
strain hardening rule, are used during the development of 
these formulations. 

In finite element technique, the equilibrium of an 
element is given by [l] 

f [B] T (dcj dv = {dF } (3.11) 

The most general egression for {da} is given by, 

id O } = [d] {d e } , for elastic case (3.12) 

and {da} = [d^ ] {d s } + {dp} for plastic case (3.13) 

since {d e } = [b] {d§} 

we have 

^ [ B ] [d] [b] dv {dr5 } = {dF } for elastic case (3.14) 

and 

* [ b ] T [D e ] [b] dv id6 } » {dF} - / [b] T {dP} dv 

v ep V 


for plastic case 


(3.15) 
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where / [B ] T [d] [b] dv or / [b] T [d ] [b ] dv is called 
V v 

element stiffness matrix and is denoted by [K]. 

In the next four sections of this chapter, during 
the development of formulations for elasto— plastic and thermo- 
plastic analysis, the main aim is to find expressions for 
[•^ep] and {dP} for different 2-D cases. 

3*4 Elgst-o-Plastic Analysis of Plane Stre ss and Axisymmetric 
Probl ems 

This analysis was first developed by Pope^ 16 ^ and 
improved by Marcal-King [33] as described in Section 3.8.2. 
Total strain is given by 

idef- = 5dS} e + *d£} p (3.16) 

where 

[_d e J = Lds x de^ d y^ J for plane stress cases 

and 

l_ds \ = |_d e ds dr de J for axisymmetric cases 

-A_ Jf X.J f 

Since 

{ dd} = [d] {d e } e 

hence using Equation (3.16) 

{do} = [d] (Cde} t - Cds} p ) (3.17,) 


Denoting the left hand side of Equation (3.6) of 
Von-Mises yield criteria by , effective stress and 


£ C < 0 x- °y> 2 + (C y-V 2 + <° e - °x » 2 + 6 T ^y] 


1/2 


(3.18) 
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for axisymmetric case. For plane stress 


Yielding begins when ‘d equals to 0 vp . 


case a = 0. 
0 


Hence increment in effective stress 


dTJ 


do + da 

3 Sc x 3 a v : 


3^ . 3^ 

d t + — — da 


31 


xy 


xy 30 @ © 


Using Equation (3.18) in Equation (3.19) and also 
usin_, che expressions for deviatoric stresses, i.e. 


s = a - c 
x x mean 


a - °Y + U Q 

x ^ 


= ? [2 


"5 I *• O “ O — O 1 

3 L x y 0 J 


and similarly 


and 


S = i [20 - O - o ] 

y 3 l y q x j 


3 © 


= I l 2o e ~ - a „l 


(■ 


We have 


“L-ls-Jtofej 


30 

" here L-wJ-Lf^ — 2 


3S^ 

2TJ 


3 T 


-2SL 


O 


3S 


’9 . 


20 


i 


and {(3oj is the vector of stress increments. 

Using Pr andt 1- Reus s flow rule i.e. Equation (3.8) 


3.19) 


U20) 

. 21 ) 

. 22 ) 


we have 
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de p V 
x Y 


\ d e p 
\ Y 


S cl 

\ 

jde p 


V 


d?v. 


xy 

P 

@ 


Calling dg P r ^ 


I s r 
I y i 

/ v 

\ T xy / 


! S 

c 


6 


(3.23) 


3 " d ?w we have Equation (3.23) as 


de 


P , 
x \ 


de p 

y 


dy p 

I xy 


de @ > 


d e p 


, s 

3 x 


\ 0 S i 
\ 1 _Z i 

I 2 * 5 ! 


V 


3 i 

a 1 

3 fe 

2 a 


hence from Equation (3.22) 


Cd e} P = de P {-££-} 

3 c J 


(3.24) 


Here d £ P is called equivalent one-dimensional plastic 


strain . 


Using Equations (3.17,), (3.21) and (3.24), we have 


= Ltv JW Cfdej* - de p f||}) 


(3.25) 


Further using one dimensional strain— hardening 
Equation (3.9) in Equation (3.25), we have 


H‘ de p = 
hence 




^30 ^ 


ds p = Lw 1 J{de J t 


(3 .26.) 


(3.26) 
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where 



By Equation (3.26) one can calculate de p , since one 
knows stress level for previous iteration and strain incre- 
ment ec^pr {ds} is calculated by finite element technique 
for c load increment . To get expression for [D ] , one 
proceeds as follows. 

Using Equations (3. 2d) and (3.24) in Equation (3.17), 

we have 

{do} = [D] ( {ds} 1 - {-|~} Lw 1 J {de} t ) 

= ^ D ep ^ {de}t (3.28) 

where 

[%] = [D]- [D]£-|fj LwJ (3.29) 

Thus the elasto-plastic analysis of plane stress and 
axisymmetric geometries, in Equation (3.13) expression for 
[p ] is given by Equation (3.29) and {dp } = {0} . 

3.5 Elasto-Plastic Analysis of Plane Strain Problems 

The formulation presented in Section 3.4, cannot be 

used for elasto-plastic analysis of plane strain problems, 

since the elimination of da i.e. the stress increment in 

z 

the third direction is to be done. Method of elimination of 
dO z was first shown in [33] , and "initial stress" solution 
procedure was suggested. In the present work same elimination 
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thod of but partial stiffness solution procedure 

h a ve be en adopt ed * 

Total strain is given by Equation (3.16) as 
t „ 


d e 


x 


Q S 


X ) 


: d 


d R.. 


dr ! 
xy j 


d e. 


'x 


de 


> + 


! r* • 

1 o 


V 


z / 


I dr 
I xy 

1 d e 


dr 


f 


xy 1 . 


d s 


( 3 ; 30) 


Jn case of plane strain problems ds t = o, i.e. 


d e e + ds p = 0 

Z 2 


In terms of stresses, this can be written 


as : 


r [ d o - v(d t 


+ d a )] + ds P = 0 

v J 


or 


da 


v (da + da ) - E d s p 

x y z 


(3.31) 


Writing first three strains of Equation (3.30), in terms of 
stresses as 


de x = E C do x - »<<% + d ° z )] + 

de y = |[ do y - 'W » 2 + do x >] + deP 


(3.32) 
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d 7 „. r = ~ dx + d-/ p 

xy G xy xy 


Solving these equations for stress increments/ in 
terms of strain increments/ we get 




1-V 


0 

h 

E 

(1 + v ) (1 - 2 v) 


1-2/ 

0 

l d v) 


0 

0 

1-2V j 
2 I 





{ 


d - v d s p — de p / 

^de t - i’ds p - de p ^ 
j Y z y 

d r fc - d^ 
xy xy 

Substituting expressions for de p , de p and de p 

x y 2 

(Equation 3.24), in the above equations, we get 


(3.33) 


f da 


x 


d e z (s + vs )de p / 

x 2c x 2 


jda 


y f 


[ p ] 


d T. 


d Sy ~ — (s + ys )de p 
y 2a y 2 


(3.34) 


xyy 


d r p - - T dl p 


xy ^ xy 

where [d] is the conventional plane strain stress- strain matrix. 

To have a similarity with the formulation presented 
in Section 3.4, let Q be an imaginary function, such that 


(3.35$ 


J § 2 . 

| = 

I - 2 -(s v + vs). 

— (S + vs ), 

_ 2 _ T j 

Lga « 


L 2 a x 2 

I'd Y 

Tf *y J 


Thus Equation (3.34) becomes 

id o} = [ D ] ( {be - d e p { §§■ } ) 


(3.36) 
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For plane strain case. Equation (3.18) becomes 

~ [<° x - o y ) 2 + (o y - o z ) 2 + (o a - cy 2 + ex^,] 1 / 2 


da . fJtaa +l!x d , + ai* d 


•g x 2 


y 7 


'xy 2 


Eliminating d o z f rom above by using Equation (3.31), we get 


do = “Z (s + DSjdo + 
2^ x z x 


(S + vs )do + 2 Y d-c 

2 O’ y 2 y V xy 


q s 2 

9 JZ 

4 U 2 


Eds 1 


{do} - | -| E di p 


( 3 . 37 ) 


Using Equations (3.9), (3.36) and (3.37), we have 


where 


l w 2 J = 


i„w 2 J{de i* 


L&JW 

+ Lf%JW IfSi+fJf 


(3.38) 


( 3 . 39 ) 


aa J 4 «2 


Substituting Equation (3.38) in. Equation (3.36), we get 
fdo} = [OgplJde } 1 

where [D ep ] • [ D] - [d] t|2 ! [wj 


(3.40) 


Hence for a plane strain elasto-plastic analysis 
expression for [D ep j in Equation (3.13) is given by Equation 
(3.40) and {dP } = {0 } . 
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3 *6 Thermo-Plastic Analysis of Plane Stress and Axisymmetric 
Problems 

Formulation presented in this section is developed • 
on similar lines, as for elasto-plastic analysis of plane 
stress and axisymmetric problems (i.e. Section 3.4), and has 
been used by several authors [12,17,19]. 

Total strain increment for a thermo-plastic structure 
is given by 

{de} t = {de f + {ds} p + { de} Th + {ds } TP (3.41) 


Th " 

where {ds} is the strain increment due to change in temp- 

TP 

erature and {de} is the strain increment due to temperature 

dependent material properties, such as E and v . 

The expression for {ds} P is again given by Equation 

Th 

(3.24) and expression for { ds} is given by 


{de} 


Th 


= dT 


<xj 

°i 

ai 


= dT {a} 


(3.42) 


where dT is the temperature increment. 


'■<CP 

Expression for {de} is obtained as follows; 


{e} = [D]" 1 {O} = [c] {<?} 


where 


[C] = 


1 

E 

E 

0 


v 

E 


1 

E 

0 


0 


1 

G 


for plane stress case 


(3.43) 
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Since 


for axisymmetric 
case 


( 3 . 44 ) 


id e y 


* gT * dT 


a f c ' 


( 3 . 45 ) 


fde} 


« v 

V - - E 

E 


- 2i + | 


for plane stress case 


{d e }' 


— £ -h — 


0 * o o; F 

V-sE V - — 
t E 


for axisymmetric case 


( 3 . 46 ) 



and 


Substituting expressions for {de} p , {de} Th 
TP 

Cde} from Equations (3.24), (3.42) and (3.45) respectively, 
in Equation (3.41), we get 

{d e } e = {de} 1 - de p {-|~-} - dT {a} - dT £ 0 } (3.47.) 


{da} 


Hence incremental stress vector {da} is given by 

= [ D ] {d e } 6 = [D] ( {d £ } fc - d e p {-§£-} - dT {a} 

a O 

“ dT ^-^{a} ) (3. 48) 


Using Equations (3.21) and (3.48), we have 

d a = b|§-J [D] ( {de} fc - ds p {|^} - dT {a} - dT {a}) 

(3.49) 

For plastically yielded point, the one dimensional 
effective stress a must be on the yield surface. Hence using 
the one dimensional strain hardening law given by Equation 
(3.9) and assuming now yield surface is a function of temper- 
ature also, we have 


da = h' de 


rP 


3 o 

Tt 


dT 


Using Equations (3.49) and (3.50), we get 
ds p = [wj ({dg} t - dT {a}- dT {d}>: 

where [w.J is given by Equation (3.27) and 


3 a 

3T 


dT 


DENW. 


DENW. 


i = H ‘ + LfrJ M { tr } 


3 0 


(3.50) 


(3.51) 


(3.52) 


Using Equation (3.52) in Equation. (3.48), we have 
expression for incremental stress vector as: 
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{do} = [D ep ] {de} t - [D ep ] (dT {a} + dT 


{O} ) 


-f 


w iHl m: 61 

DENW.. 


(3.53) 


where [d ] is given by Equation (3.29). 

Cornparing Equation (3.53) with general expression 
for {do} given by Equation (3.13) , one gets expression for 
f D ep] as given by Equation (3.29), and expression for {<r).P } as 


Cd.p} 



■(dT {G.}- + dT 


9 f~C~ 
3T 


io } ) 


+ 


[D] f 


9 0 3 C 
3 O 3T 
DENW 


dT 


(3.54) 

This can be noted that for thermo-plastic analysis, 
further load increment depends upon stress field from the 
previous iteration and hence it is absolutely necessary to 
calculate stress field very accurately. 

3 . 7/ Thermo- Plastic Analysis of Plane Strain Problems 


The welding of long plates or cylinders is essenti- 
ally a plane strain case, as explained in Chapter 6. In [ 20] 
this case has been simplified by postulating an expression 
for de" 13 . Some authors, such as in [7,17;], keep a separate 
identity of stress increment in z-direction for this case 
also. In present work, the elimination concept of Section 
3.5 for plane strain elasto-plastic analysis has been exten- 
ded for thermo-plastic analysis of plane strain problems. 

This again makes stiffness as well as further load increment 
a function of previous stress level and hence determination 
of stress level in every step very accurately is the must. 
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as small error in stress level keeps propagating and diverges 
through further iterations. 

The total strain increment . 

is given by Equation (3.41), which in e;xpanded form 
is 


d c 


de 


t 


< . 


a. r ! 

XT’ -j 

j 

[ d6 z I 


de 


x 


\ 


ds. 


y 

xy j 


< - 1. 

1 dr 


d e 


z J 


fd e n 
x 1 


ids 


d y 


xy 


d e. 


J 


+ 


For a plane strain case de u 

z 


de 


x 


Th 


d s 


y 


I dy 


xy | 


d e 


i 

J 


- 0, hence 


de^T + de P + de Th -f de TP = 0 

z z z z 


i d e 
i x 


d s 


+ 


dy 


xy 


d e 


(3.55) 


(3.56) 


Substituting de in terms of stresses from Hooke's 

Z 


law, we have 


1 


v l-o " ^ (do + do ) 1 + de v + de + de 
E u z x y J z z z 

or do 


= 0 


v (d a + do ) - E(de p + de Th + de TP ) 
x y z z z 


(3.57) 


Also rest of the relations of Equation (3.55) can be written 
in terras of stresses as 


de 


t 

x 


f d a - v (do + do ) 1 + de K + de Th + de r 
L x y z J x x x 


TP 


E 


de v ' = i [do - y (do + do ..) 1 + d e p + de Th + de 
y E L y z X J y y y 


TP 


(3.58) 


d y 


3cy 


i dr + dr p + dy TP 
G xy xy xy 


(*.* dr Th = 0) 

xy 


TP 
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Substituting for dc z from Equation (3.57) and solv- 
ing foi incremental stresses in terms of incremental strains , 
we get 


■d a 


x 


<da > 

? v 


!d 


(1 + ^ ) (i -'2 v) 


de 


dr 


x 

t 

y ’ 

t 

x< 


X 

(de p 


1-y 

V 

0 

V 

1-v 

0 

0 

0 

l-2r> 

2 


_Th 
as - 

X 

de TP j 
x ! 

!) - 

ds Th _ 

y 

. tp ! 

de : 

y > 


dr p - d r TP 


xy 


xy 


J 


|1 

- > ji > 

k 

Substituting expressions for de p ds p and de p 

x y z 

from Equation (3.24; , Equation (3.59) becomes 


(3.59) 


id a 


x 


/da 


'xy 


i 

de 1 " - 

_ 3_ 

(s v 

+ V 

s)dl p - 

de Th - 

de TP j 

1 

X 

2 a 



Z ’ 

X 

X 

•r^~ 

ii 

r — i 

M 

de fc - 

__3_ 

(s 

+ 

S )de p - 

de Th - 

de TP , 

/ \ 

y 

2 a 

y 


2 

y 

y 


dr fc - 

3 

X 

de p 

- dr TP 




xy 

a 

xy 


xy 



i 







.... ^ 


Ev . , Th , TP X 

(1+2?) (1 - 2 v ) (Qe z + ds z ) 


(3.60) 


where 
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$ j 1 

i A } = j 1 > 

I 

|° 

L J* 

Using the notation defined in Equation (3. 35)/ we have 

* da ! = [D] (fde) 1 - dsP{|§}- ide} Th -!de] TP ) 

EV , „ Th TP 

“ TTT V)(l- 2V) (de z + d£ z ) * A J (3.61) 

Using Equations (3.50), (3.38) and (3.61) , we have 


e p = |_ W 2 J < * d e ~ ids } Th - { 


d s} TP ) 


(■ 


EP 


(1 + V ) (1 - 2 V ) 1-3 


^ G 


DENW, 


30 

aX 


dT 


DENW. 


(de™ + de TP, 

Z Z 


(3.62) 


where |w_J is defined by Equation (3.40) and 


DENW 2 = H' + ijflj [D] !§*} + f 


f 3Q. 


9 z 


(3.63) 


Using Equation (3.62) in Equation (3.61),. we have 
Cd°! = [%] (Cde; 11 - !de} Th -f dej TP ) + Xl i fe?. ^ — 


DENW, 


irfvrU w ISJ 


3 fz 
3Q , 2 a 


v 


DENW„ 


ii-iv) f A J] < de 2 h+ ds r> 


, TP. 
Z 


(3.64) 

Again comparing Equations (3.13) and (3.64), we have 
expression for [ D ep ] as given by Equation (3.15) and 
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{dP} = » [D ep ] ([1 1 0J T a dT + ^JL^J dT { c } ) 




- s 

3 z 


E 

(1 +V) 


L- Co] t&} 


V 


3 c DENW„ 


DENW„ 


(i - 2V ) fA3]Cdef + de^; 

(3.65) 


121 


TP, 


3 ' 8 H iterative Solution Algorithm 

Four sets of formulations presented in Sections 3.4 
through 3.7 can be used for elasto-plastic and thermo-plastic 
analysis of 2-D structures. Different solution procedures 
are available in literatures. An iterative solution algorithm 
adopt'd in this work, for solving nonlinear equations, is 
presented here [6,27,33]. 

3.8.1 A Solution Step 

A typical step of the solution is described with the 
aid of .Figure 3.5(a), which is a single degree of freedom 
representation of an actual multi degree of freedom problem. 
Let under loads {R } A , the correct displacements 16 } A and 
structure tangent stiffness [kJ a are known. Displacements 
produced by the next load increment are computed as 



Ma 

* a<5 -*ab * ar -*ab ~ * r *b “ * r *a 

(3.66) 

and 

£sj b 

= £6} A + 

(3.67) 


Stresses, strains and tangent stiffness must now be 
updated. Displacements {Ad} of an element are extracted 
from £46}^* Then, for each sampling point in the element. 
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total strain increment {As} =[b] {Ad} is computed. Once 
{As} are known following quantities are computed from diff- 
erent expressions cerived earlier. 

Expressions for ds^ for all the four cases considered 
in Sections 3.4 to 3.7 can be written in the form 


de P = [f x [ {d e} t + f dT 


(3.68) 


where for plane stress and axisymmetric elasto— plastic case 

L f lJ = L w ii and f 2 = 0 (3.69) 


for plane strain , elasto-piastic case 

L f i J = Lw J and 


f 2 = 0 


for plane stress and axisymmetric thermo-plastic case 

L f lJ = L W lJ 


(3.70) 


and 


f 2 " 


- L w iJ <{«} + ?a}) - 


do 
d t 

DENW, 


for plane strain thermo-plastic case 


[f 1 J = l,W. | and 


(3.71) 


' L w aJ ( £°i + a a£ L WJ) 


( 


EV 


Cl + v) (1 - '2 V ) U C j {A 5 + 2 a 


DENW, 


( a — i (V - ~ E)(G + O ) - C ^ ) 
^ E v E x y „2 z ; 


3 Q 
3 T 


DENW, 


(3.72) 


Using Equation (3.6a), we have 
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AeP 


b 

' L £ J < de } t + 


b 

/ 


f 2 dT 



[As } ^ 


f 2a A ' T 


(3.73) 


Also from the flow rule given by Equation (3.24) 


fA V 


3 0 

f {- } dgP co ? ! g 7 , - 

a ac - *~ 3 a As p 


(3.74) 


Also expression for {d<r} for aJ 1 +-Ho -p 

ci ror an the four cases consi- 
dered in Sections 3.4 to 3 i, r>a n Ho 

° J * 7 ' can be written in the form 


{da } _ [ Dep ] {de } b + {f } dT 


(3.75) 


where for plane stress, plane strain and axisymmetric elasto- 
plastic cases 


ff 3 } = £0J 


for plane stress and axisymmetric thermo-plastic cases 
if 3 J = * [ D ep] ( £ a ) + la} ) + ^LilZLlZ 


(3.76) 


for plane strain thermo-plastic case 


ff 3l = - [%] < !aj + [OJ , + 


(1 + V) C" L D ] fft) 


c»] 


o S 
3 z 


denw. 


3Q 90 
DENW_ 


0> 


(3.77) 


'3 0 J EEiW 2 71 _ 23 ?) {A * 

o • 

f y _ — TTWcr 

- E 

So using Equation (3.7,5), we have 


+ C ) - E- 0 2, 
■ E z 


(3.78) 
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b 

{Ac}= / [D ep ] {de} t + / { f 3 } dT 

^- D ep-i a ^ + a dT (3.79) 

Then values of £ p/ { e p } and {o} are updated by adding 
corresponding increment to their previous values. The up- 
dated c can be computed from the updated {c} . Also from 
the updated e p/ H' can be calculated from the input material 
properties. Now for the current stress level and material 
properties, [ D e p] snb £dP} are calculated from the different 
equations given above depending upon the case. Hence updated 
stiffness matrix and load increment for the point B in 

3.5(a) are calculated. In this solution - algorithm 
one has to keep records of e p , U p }, c and {a} at each 
sampling point in each element and store the updated { 5 } if 
structure displacements are needed. 


3.8.2 Elastic-Plastic Transition 

When a sampling point makes the transition from 
elastic to plastic within a single load increment, one has 
to switch on from equations of elasticity to equations of 
plasticity. Hence a correction is needed. In the present 
work the method adopted is as suggested in [36] . In this 
method when a sampling point goes to plastic state in a 
single load increment, a factor *m? is calculated by 




(3.80) 
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where o is the previous elastic stress and a is the 
OXQ new 

present elastic stress which is more than yield stress a . 

X hr 

Then [D ] is calculated by 

[ D ep ] = [D] - (1 - m) [ D] {-||j |_W X J (3.81) 

for plane stress and axisymmetric cases and 

[%3 = M - (1 - rn> [D] £||} l_„ 2 j (3.82) 

for plane strain case. 

Hence (1-m) represents a fraction of partial yielding 

of sampling point. This [d 1 is used to calculate stiffness 

matrix. Displacements calculated by this [D ], will be 

nearer to the yield surface, say point 1 C' in Figure 3.5(b). 

Hence more iterations may be required to bring all the 

stresses on the yield surface, by perturbing value of 1 m‘ 

accordingly. After the convergence. Equations (3.68) to (3.79) 

are used to calculate different parameters, replacing 
+- Th Til 

by (1-m) As and also As by (1-m) As . It may be noted 
that this replacement is different than suggested in Cook [3], 
where replacement Ae by (1-m) As is suggested only. Present 
author found this replacement is incorrect, since pseudo 
thermal strain should also be subtracted from total strain, 
before multiplying by the fraction of partial yielding. Also 
while calculating unbalance plastic load {dPj , this vector 
should be again multiplied by (1-m) . 
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3.8.3 Calculation of Unbalance Forces 

A drift appears, when upper solution procedure is 
adopted, as shown in Figure 3.5(a). To avoid this, following- 
method has been used. 

Upon reaching R fi , the computed displacement < 5 Q is 
smaller than the true value, hence the resulting stresses 
are smaller too and are not adequate to resist the applied 
load Rg . If a corrective load AR c were added to AR, the 
computed 6 B would be larger and more nearly the true value. 

The corrective load, which is also called unbalanced force, 
is the difference between applied loads and resistance pro- 
duced corresponding to the calculated displacements. Hence 

= fR}- E / [B] T {c} dV (3.83) 

COr e Vol 

where the summation extends overall elements of the struc- 
ture. Forces £R} and stresses {a} are those at the end of 
a step. Correction is applied in the subsequent 

step and hence equilibrium equations for step B to C is 

Mb * aa i BC = {ar! bc + !AR! cor,B 

There are other methods available such as Newton 
Raphson and modified Newton Raphson, to avoid drifts. These 
methods ask for more number of internal iterations, within a 
load increment. But in the present work it is seen that 
just by adding unbalance force to the next iteration, it is 
possible to follow yield surface closely, provided load 
increment is kept below a particular value. For example, 
in case of thermal load, temperature increment of a sampling 
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should not be more than 50° C in one load increment 
step. But it is necessary to do internal iterations to bring 
a sampling point on the yield surface, whenever there is a 
transition from elastic to plastic stage, as explained in 
Section 3.6.2. 

3 * 9 A Computer Program and Chech Out Problems 

Depending upon the formulation presented in Section 
3.4, present author developed a computer program PLSTIC 
before this work and was reported along with the check out 
problems in [27] . During the course of this work same 
program was modified for elastic-plastic plane strain case 
(formulation 3.5), thermo-plastic plane stress and axisymmetric 
case (formulation 3,6) and thermo-plastic plane strain case 
(formulation 3.7.). This computer program THESIS (THErmo- 
plastic analy sis ) was tested for number of sample problems to 
check its capability. Some of these sample problems are 
presented here to demonstrate some of the 

analysis capabilities of the program THESIS and to show the 
accuracy of the solutions and hence comparisons are made with 
the results of other authors, wherever available. 

3.9.1 Thick Cylinder Under Internal Pressure ( A Plane 
Strain Case ) 

This problem was solved to test elasto-plastic plane 
strain analysis capability (corresponds to formulation 3.5) 
of THESIS. For this a cylinder under internal pressure was 
solved for plane strain elasto-plastic case. The cylinder 
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FIG- 3-7 RADIAL STRESS VARIATION OF E LASTO-PL ASTiC 
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was characterized by the following non-dimensional parameters: 

outer radius/inner radius = 2 

v = 0.3 

yield stress/E = 0.002 

In Figure 3.6 finite element discretization is shown 
for this geometry and Figures 3.7 to 3.9 show the variations c 
°q sna with radius, when plastic zone radius is 
equal to 1.5 times the inner radius. In the same figures, 
a comparison has been made with the results quoted for the 
same problem in "[8]. Figures show a good agreement of differ- 
ent parameters except near to the inner radius, where present 
program computes a much better results than in [8]. This can 
judged from the variation of a r , where a r must be equal to 
-p at the inner surface to satisfy natural boundary condition. 
Ref. [8 ] failed to reproduce this boundary condition, whereas 
present program followed this condition with smooth variations 

of cj , o Q and 0 . 
r d z 


3.9.2 Thick Cylinder Under Radial Temperature Gradient 
( Plane Stress as well as Axisymmetric Case ) 

The results for this problem are presented here to 
show the solving capability of program THESIS for thermo- 
plastic analysis of plane stress and axisymmetric bodies. 

Theoretical solution of this problem is quoted in 
[2]. The problem is characterised as a thick cylinder of 
2.25 inner to outer radii ratio, is under radial temperature 

gradient. The radial temperature drop is increased from 0 °C 

E a A T 


to a value for which 
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Since essentia11 * this problem is a one dioensional 

case, hence it has been solved by thesis as a plane stress 
as well as a*i synmetrio oase . ^ ^ ^ 

the plane stress and asymmetric finite element discretize- 

tion of this problem respectivelv 

pectively. Results obtained from 

both the analysis are similar. 

in Figure 3.12, progress of plastic front is shown 

as a function of upto the value of - 6 0 f 

YP ~ - o.o, for 

which results are quoted in f 2 1 same ^ • 

L i' bame problem is solved in 

[12 ] and Figure 3 . 12 also shows results from [121 it is 

.... M J r 12 ^ 

With theoretical result, than quoted in [ 12 ]. In Figure 3.13. 
variations of O r and along radius are shown for gg AT = 

4.525 and a comparison has also been made with theoretical 
results from [2]. 

3 ' 9 ' 3 cradw 

Lastly to check the thermo-plastic analysis capabi- 
lity for plane strain problems of THESIS , solution of thick 
cylinder under radial temperature gradient , considered as 
plane strain case is presented here. Same problem is solved 
m [ 35 ] and those results are also shown for a comparisons. 

The problem is characterised as follows: 
inner radius/outer radius = 4 

E = 31 


v 


0.25 



AT Q Eor 

°yp(1-V) 
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FIG-3-14 THE GROWTH OF PLASTIC ZONES IN A PLANE 
STRAIN CYLINDER UNDER RADIAL TEMPERAT- 
URE GRADIENT 
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Progress of plastic front is shown in Figure 
3.14. A comparison with [35] shows both results are not in 
good agreement. Present author feels this is essentially 
because, yield criteria used in [35] is Tresca and in THESIS 
yield criteria is Von— Rises. Similar differences in progress 
of plastic front because of different yield criteria are 
reported in [6] for only elasto-plastic analysis. 

In Figure 3.15, variation of 0^ along radius is 
shown, when plastic front progressed upto 1.5 times inner 
radius . 



CHAPTER 4 


™l . ^°~!^ . gHANlCA L ANALYSIS OF 2-D STRUCTURES WITH CREEP 


4.1 Introduction 


4.2 Creep Laws and Representative Incremental Creep Strain 
Form ^ 


4.3 Derivations of [c ] and fds c } atrices for Different 
Creep Laws 


4-3.1 Time-Hardening Creep Law 


4.3.2 Power-Creep Law 


4.3.3 Exponential-Creep Law 


4.3.4 Strain-Hardening-Creep Law 


4.4 Thermo-Mechanical Analysis of Plane Strain Case with 
Creep 


4.5 Solved Examples 


4.5.1 A Note on Material Properties 


4.5.2 One- Dimensional Creep Stress Relaxation at 

Constant Temperature 


4.5.3 One-Dimensional Creep Stress Relaxation* for 
Varying Temperature 


4.5.4 One- Dimensional Creep Stress Relaxation for 

a Plane Strain Case 



4 • 1 Introduction 


In order to investigate stress annealing behaviors 

creep 

o we ed structures, the effect of /strain should be included 
in the theory presented in Chapter 3 . To consider creep 
strain Cyr and Tete^ 28 ^ suggested separation of creep strain 
from elastic and plastic strains in the formulation. In [is] 
authors adopted this way of including creep strain in exis- 
ting elastic-plastic formulations of tangent stiffness method. 
They derived incremental stress-strain relations for several 
creep laws, such as power, exponential, time-hardening and 

h ar< 3 en i n< ~f • These incremental stress— strain relations 
are used here for analysing annealing of welded structures. 

These four creep laws are given in Section 4.2. 
Derivations of incremental strain, as done in [18J, are 
presented in' Section 4.3. A modification for plane strain 
case including creep strain, as done for thermo-mechanical 
analysis in Chapter 3, has been done by the present author 
and is presented in Section 4.4. In Section 4.5 solutions 
of some sample problems are shown for verification of this 
formulation. 

4.2 Creep Laws and Representative Incremental Creep Strain 
Form 

Many mathematical relations are used to represent 
creep behavior of a material. But four creep laws, such as 
power, exponential, time— hardening and strain— hardening seem 
to be popular and a combination of these laws can represent 
the creep behavior of a material for entire range of annea- 
ling temperature of a welded structure successfully. Hence 
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all these creep laws are used in the formulation adopted in 
this work and any of these laws can be used for any temper- 
ature range depending upon the user's choice. 

In one— dimensional stress state, these four laws can 
be represented mathematically as 


6 C = 

J3o n 

* 

• 

power-creep law 

(4.1) 

e c = 

*c 

A exp (B 0 ) 

„ _ y , m-1 


exponential-creep law 

(4.2) 

E = 

m a 0 t 

* 

time-hardening-creep law 

(4.3) 

e c » 

ra a 1/m ' o r /m 

1- - 

(e C ) m : 

strain-hardening-creep 

law 

(4.4) 

where 

o = axial stress; 

strain rate i.e. 

c 

e = creep 

de c 

dt 

> axial strain; £ C = creep 



p , n, A, B, m, a, r are all constants and functions of mate- 
rial properties. 

In three-dimensional stress state, assuming that the. 
creep behavior is described by the theory of incremental 
strain of plasticity , the component of the creep strain 
rate is proportional to deviatoric stress. Hence 

ie C ] = — ' {4 ‘ 5) 

20 

where f C : proportionality factor, called equivalent 
creep strain rate 

a : equivalent stress as defined earlier 

L S J = L s xWzJ 
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Taking account of the effect of variation in the 
state of stress and creep properties , the general form of 

creep strain increment {de} cr during a time increment dt 
is given by 

{ds} cr = [C c ] { do} + {de G } (4.6] 

Hence for different creep laws matrices f C ] and 
{d s } should be derived in terms of state of stress and 
creep properties. 

^ ^ D erivations of [c„1 and {d£ G } Matrices for Different 
Creep Laws ~ 

Comparing Equations (4.1) and (4.3), we find that 
if 

r = n ? a = £ and m = 1 (4.7) 

then Equation (4.3) reduces to (4.1). Hence it can be said 
that power-creep law is a particular case of time-hardening 
creep law. So form of [C c J and {de c } for time-hardening 
creep law can be used by substituting equivalence given in 
Equation (4.7) for power-creep law. 

4.3.1 Time-Hardening Creep Law 

Substituting for 1-D creep strain rate in Equation 
(4.5) from Equation (4.3), we have 

ii C } = | m a o 7 ' 1 t™" 1 {S} (4.8) 

Using the notations 

b = m a and M = m - 1 (4.9) 



We have 
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Assuming that b, 5 and {S} in Equation (4.10) vary 
linearly and Jl , v are constants over the time interval dt , 
and integrating this creep strain rate during the time inter- 

val fc l i i + dt ' wh ere t ][ is real time, the vector {d£ c } 
and matrix [C c ] can be expressed as follows 


Cde c } = ^ (K^b + db) {s} 

and 

[C c ] = ( Kl b + k 2 ab)c r - 1 [c ±J ] 

where 

K 1 = * (t i + dt) M+1 - ( tl ) M+1 }/(ju+ 1) 

K 2 = (t l + dt )^ + V(M + 1) - {(t ± + dt) M+2 
- (t 1 ) 4+2 }/{(,U + l)(ju + 2) dt} 

K 3 = (t 1 + dt) ,1+1 /(}l + 1) - 2(t x + dt) M+2 / 

£ (M + 1 ) ( JUL + 2 )dt } + 2 { (t x + dt f' +3 - 
(t 1 ) M+3 }/{( n+ 1)(UL + 2) (11 + 3 ) ( dt ) 2 } 

Also 


11 

== 

F S. S + 
X X 

1 

C 12 

= F S S - 0.5 
x y 

13 

$s 

2F S T 

x xy 


C 14 

= F S S - 0.5 
X z 

22 

s= 

P S y S y + 

1 

C 23 

= 2F Sy 

'24 

= 

F S y S 2 * 

0.5 

C 33 

+ 

P 

S’ 

P4 

KTf* 

II 

34 

= 

2F s z T xy 


C 44 

= F s z S 2 + 1 


(4.11) 


(4.12) 


(4.13) 


(4.14) 


(4.15) 


(4.16) 
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where 


ij 


- C (i/j = 1 ~ 4) and f = 9 (r- i)/(4c 2 ) 


4.3.2 Power-Creep Law 

It is apparent that for the power creep law, the 
creep strain increment is obtained by putting b, y and ji in 
Equations (4.11) through (4.16) equal to P , n and zero, 
respectively . 


4.3.3 Exponential-Creep Law 

Substituting for one-dimensional creep strain rate 
from Equation (4.2) into Equation (4.5), we have 


r * G 

t 8 i 


2 c 


A exp (BO ) {S} 


(4.17) 


Again assuming variations of A, a and {S} as linear 
and no variation of B, over the interval dt, we have 


{d e c } = — (A + ~) exp (bo ) {S} 

20 


[c e ] - * <f ♦ Ccij] 


(4.18) 


(4.19) 


Such that 


C 

C 

C 


1 1 

11 

= G S S + 
XX 

l 

C h 

= G S S - 0.5 
x y 

> 1 

13 

= 2G S x T xy 


C i4 

= G S x S z - 0.5 

22 

“ G S Y s y + 

l 

c 

u 23 

= 20 S y T xy 

«t 

’24 

= G S S_ - 

y 2 

0.5 

c* 

C 33 

= 4G T T 4* : 

xy xy 

'34 

= 2G S„ T 
z 


C’ 

C 44 

= G S S 4-1 
z z 


(4.20) 



also 


C ij ~ C ji j = 1 - 4) ; 
G = 9 (BO - 1 )/(4 *5 2 ) 


and dA is the change in A due to temperature increment in 
time dt. 

When Equations (4.18) to (4.20) are used for analysis, 
a time, increment dt should be chosen to satisfy the ineq uali ty 
B do < 1, as an exponential function is approximated by a 

series of polynomials in the process of derivations of these 
expressions . 


Strain-Hardening Creep Law 

In case of strain-hardening creep law, given by 

Equation (4.4), material constant m is generally less than 

unity . Since at the beginning of the iterations creep strain 
c . 

s is zero, hence from Equation (4.4), creep strain rate e c 
is infinity. So, direct application of this equation is not 
possible for computation. A method suggested in [18], has 
been used in the present work for this creep law implemen- 
tation . 

In [18], a new creep law, called total-strain creep 
law, has been defined and is given by 

s G = a a r t m (4.21) 

Hence strain-hardening creep law (Equation 4.4) is 
now a solution of simultaneous equations of Equations (4.3) 
and (4.21) for parameter t. In such a case where the 
stresses in a continuum do not change in a specified period 
of time (t^ < t < t 2 ) , the following method is completely 



equivalent to the direct- anr-ii-s,- =»*.-! r 

rect application of Equation (4.4) to 

evaluate the creep strain. The creep strains are calculated 

by integrating Equation (4.3) for the specified interval 

^ fc eq ^ ~ ^eq + ^2 “ 4n which the equivalent time 

t eq 13 calculated by Equation (4.21) of a function of stres- 
ses and creep strains at the real time t^ in actual stress 
relief annealing, the stresses in a structure vary during 
the treatment. The changes of the stresses during a small 
time interval can be regarded as invariant. In this case, 
the creep strain increment under the strain-hardening creep 
l3.w can b<3 evaluated by the above-mentioned method. 

In the three-dimensional stress state, the equivalent 
time t e g is given by the equation 

t eq = U C /(af)} 1/m ' (4.22) 

Hence Equation (4.3) must be integrated for analysis 

during a time increment dt(t__ <t<t_+dt). So, the 

eq eq 

results given by Equations (4.11) through (4.16.) can be used 
by replacing t^ by t . Also at the first incremental step, 
t is equal to zero because there is no creep strain* 

tsq 

4-4 Thermo- Mechanical Analysis of Plane Strain Case with 

Creep 

In case of elasto-plastic and thermo-plastic anal- 
ysis of plane-stress and axisymmetric structures with creep, 
formulations presented in Sections (3.2) and (3.4) can be 
used, after adding one more strain terms corresponding to 
creep, given by Equation (4.6). Hence no attempt has been 
made to rederive those formulations. But in case of plane 
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strain structures, an attempt being made here to rederive 
the formulations presented in (3.3) and (3.5), as many new 
vectors appear in the final expressions if creep is consid- 
ered. Derivation given below is essentially for thermo- 
mechanical analysis, which of course is a generalisation of 
elasto-plastic case. This conversion can be done by dropping. 
The components of {d^}, {de} Th and {de} TP from the final 
expressions . 

lotal strain increment is given by (considering creep 

also) 

■ t* 0 Ti r pVv T* P 

{de} = [de] + {de} p + {dej + {de} + {ds} cr (4.23) 

Substituting |d s} G in terms of stress increments 

* 

from Hooke's law and {de} cr from Equation (4.6) in Equation 

(4.23), we have 



1 

— V 

0 




1 

' I 

** V 

1 

0 

-V 

{da} + {de} p + 


0 

0 

E 

G 

0 

- 

- 





+ {d S} TP + [C c ] 

id a} 


-y 

-V 

0 

1 

- 



** 




+ {ds c 3 

(4.24) 

Denoting 

components 

of 

matrix [ C c ] as c n' 

C 12 etc " 


we rewrite Equation (4.24) as 
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fde }' 


! +C U 


■"^2 c i 


E^l 1^22 


-V 

E ^14 


-V 

E ^24 


~ +C c 

0 33 c 34 


E-Sl =£^2 C 43 i+c^ 


[doi + ide}P tide)™ +{deJ TP + (d $ c ) 


(4.25) 


Since for plane strain case de^ = 0, hence from last equation 

of Equation (4. 25), we have 


mm'P 

( “E + c 41 )d0 x + {Z E + C A0 )da ,r + + (h + C„ JdO 


tr* 1 vja/ww v- m, -r V7T t U „ . , 

k 42 y 43 xy E 44 


+ de^ + de Th + de TP + ds c = 0 

z z z z 


(4.26) 


T T Te C ■ , ) 5^- E c 41 )d0 x + (1 " E C 42 )d0 y- 


E C 43 dT ™ 

43 xy 


E(d©P + de^ h + dz TP + dg®)} (4.27) 


Hence eliminating d o- from first three equations of 

z 

Equation (4.25) and also using flow rule 


{ d e}* 


3 dj> 
2 a 


(4.28) 
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We have 


{d zT 


[c] {d0 } + de p {|2, + lde} Th + j dsJ TP 


+ {de G } 


where 


Th 


+ CH 3 (det 11 + de TP + de G ) 

Z z Z 


(4.29) 




1+EC 


(* 


11 


E 

V +EC 


)+( 


W-EC 


41, 


1+EC 


44 


14 


E 


Symmetric 


-W+EC V-EC . 

( — ~ — ~)+( ^2 


E 

-v+Ec 


l+EC 


44 


)-- c -( 


'43- 


13 1+EC 


)* 


44 


(• 


14, 


E 


(-W+EC 14 ) 


1+EC 00 V-EC,_ C 

(_-22) + (__i2)-x- c (r _43 


E 

V+EC 


'1+EC 


44 


'23 1+EC 




44 


24 


E 


) 




” +C_ «*”C_ . ’ 
G 33 34 


EC 


43 


1+EC 44 


(4.30) 


{ 3Q } 

l 30 i 


3 V “ E C 14 

£— (S + -•■■ ■ — s ) 

2 1x1 +E C 44 2 ’ 


3. (s + * ~ E 2.24 s } , 

*■ (S y + 1 +E C J4 z < 


20 


“ (T 

a W 


'44 


E C 


34 


(1 + B C 44 ) 2 


(4.31) 


and 
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{H} 



(4.32) 


,i-l 


Denoting [ Cj = [D], we have from Equation (4.29) 

{da } . [D] ({dej* - d«P{|5j- {d e} Th - { de} TP - { dl c 3 ) 


a 

Hence 

da = 


- [5 ]£HJ (def +de TP +d %c ) 

2 Z 2 

Using Von-Mises yield criteria,, we have 
[(a - a) 2 + (a - a) 2 + (a, - a,) 2 ^ 


(4.33) 


sf 2 " x y ‘ ’z “x' ’ - y w z ‘ - X y. 


2 , 1/2 


6 x ] 
XV J 


s 

- — do + ~ ^ da + — 2 


S, 

a “• 'd 

Using expression for da from Equation (4.27) in 


2 x ‘ 2 v? y z ~ xy 


3S_ 3t 

da + — dt 
2a 2 a 


(4.34) 


(4.34) we have (also noting Equation (4.31) and Equation (3.50)) 

da = h’ de p + ^ dT = i ~ | fd a 3 - ~ -2. — : ■ — - 

3T LaaJ 1 ^ 2 *5 (1 + E C 44 ) 


(de p + de Th + de TP + d£ c ) 

v Z Z 2 Z 


(4.35) 


Solving for dS P from Equation (4.35), we have 

de p « [w J (Cde} t - Cde } Th - {de } TP - {de c } ) 
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E 


9 a ' J T 2 ~ 7 t ~F 7 1 — r } 

z c u + E C, .) J 


a g 

3T 


dT 


DENW 


DENW 


— <d E f + de TP + d|C 

z Z 2 


(4.36.) 


where 


L|f.| C 5 ] < Wef - ( de} Th_ fdeJ TP_ (dS c r) _ 


9Q 

a o 


J [D] {H ] (d^ h + de f + d |°)- 1 J 


7 S 

3 2 


2 ^ U + E C, .,) 
-±*± 


Wj: 


and 


(de 


Th x . TP . a C , 3 “ 

2 + ds 2 + d e? ) - 1^- dT} 


DENW 


(4.37) 


DENW 


9Q i r~i f 9Q ' . 9 z 


H ' + l5f J[o]tf^J + 


30 4 ~2 (1 + E C 44 ) 


(4 • 38 ) 


Substituting expression for di p from Equation (4.36) 

into Equation (4.33), we have 


{da } 


[D ep ] ({de}*- {de} Th - { de } TP - { d &- c } ) + 


!£. 


dT 


[D] { |2}^Il + [Dj({f|}LiJ{H} 


DENW 


3 s * 


E 


+ |f} .■ £ . - {Hi ) (de^ h + def + d e ^) 


DENW 


(4.39$ 


where 


Hence in Equation (3.13), „e have [ D ^] s [5 ] and 


»P} = [D ep ] (!de} Th +{ del TP + td *c }) 


n a 

DT 1 


DENW 


[ 5 ]'|§>-[E] (iff I 1_WJ [H 


_ s 

3 2 


E 


Apr 2 1 + E c , , 

+ J ^ r i4 

DENW 


~ {H } ) (de Th + de J " t ' + de c ) 
z z z 


TP 


(4.41) 


It can be seen here explicit egression for [D] 
is very difficult to find, and hence one has to compute matrix 
[cj from Equation (4.30) and then invert it numerically for 
all gauss points for every iteration. This inversion is 
likely to be easier from numerical point of view, as [C] matrix 
will be a diagonally dominated matrix. 

4 . 5 Solved Examples 


Program THESIS described in Chapter 3, is then modi- 
fied for taking into consideration creep effect, depending 
upon the formulation presented in Section (4.4) . This program 
is then tested for different sample problems.' Some of these 
are presented here, also to show the accuracy a comparison 
is made with the results of same problems from [18] • 


4.5.1 A Note on Material Properties 

Material properties used to solve different sample 
problems are same as given in [l8], which are essentially 
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for quenched and tempered high strength steel H.T. 60. These 
properties are shown in Figures 4.1 and 4.2. This material 
seems to obey power creep law above 550°C and exponential 
creep law below 550°C. 

4.5.2 Qne~ Dimensional Creep Stress Relax ation at Constant 

Temperature ’ 

The problem consists of a 1-D bar, carrying axial 
initial stress of 12 Kg/mm 2 . Bar is fixed from both the 
sides and kept at a constant temperature of 600° c, for 3 hours 
Since temperature of the bar is above 550°C, hence 
as per the material properties described here, it follows 
power creep law. Using the properties as 

n « 4.2 and j3 - 2.7 x 10*" 9 

Since program THESIS is essentially for 2-D struc- 
tures, hence present problem is treated as plane stress case. 

The results are 

shown in Figure 4.3, A comparison has been made with results 
shown in [ 18 ] and found in excellent agreement. 

4.5*3 One- Dimensional Creep Stress Relaxat i on for Varyin g 

Temperatures 

In conventional stress relief annealing, it is a 
common process that the structure is gradually heated up to 
a holding temperature. So, the initial stress decreases 
during the heating stage by reduction of Young' s modulus and 
the yield stress, and also due to increase in creep strain. 



ESS VARIATIONS DURING HEATING AND C 
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To study this / problem defined in Section 4.5.2, is 
solved for a heating rate of 600° C/hr for one hour and then 
holding at this temperature for 2 hours and then cooled down 
to room temperature at 600° c/hr. 

'Ihe material properties used are as defined in Section 
4,5.1. Ihe x esults are shown in Figure 4.4 and also compared 
for first 3 hours with results quoted in [l3j and found in good 
agreement. Results show at the range below 400°C, stress decr- 
eases due to decrease in Young 1 s modulus and at the range above 
400° C the stress indicates lower value than one estimates from 
only the change of Young's modulus. In this case, the stress 
is in the elastic range even if high temperature is held. Then, 
it decreases not only due to decrease in yield stress but also 
because of generation of creep strain. 

4 .5.4 One-Dimensional Creep Stress Relaxation for a Plane 
Strain Case ' 

Example solved in Section 4.5.3, is also solved as a 

2 

plane strain case for initial ^ =12 Kg/mm and no inplane 
stresses. Though results of this problem is not available in 
literature, but author ejects to obtain a similar variation of 
stresses as was in Section 4.5.3, since problem essentially is 
1-D. Also since while solving as a plane stress case, bar is 
assumed to be of small and constant length, and hence also can 
be solved as a plane strain case, assuming strain is zero along 
the length of the bar. Results of this example are shown in 
Figure 4.5, and found as per the expectation. 
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chapter 5 

THERMAL RESPONSE ANALYSIS FOR NONLINEAR 
HEAT TRANSFER PROBLEM 

5 . 1 Introduction 

5.2 Theory of Two Dimensional Linear Heat Transfer 

5.3 Solutions of Differential Equations for Transient 

Analysis 

5.3.1 Finite Dif f erence 

5.4 Nonlinearity in Heat Transfer 


5.4.1 Change of Material Properties with 
Temperatures 

5.4.2 Consideration of Thermal Radiation 


5.4.3 Phase Change 
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5.1 Introduction 

Calculation of residual stresses during different 
mechanical processes/ such as welding , casting/ grinding 
etc./ necessitates prediction of temperature distribution 
with a high degree of accuracy . Since present work is aimed 
to determine residual stresses during mechanical processes/ 
a computer program is required to calculate temperature dis- 
tribution for nonlinear heat transfer problems. A computer 
program WELTEM was developed by the present author, as a 
part of this integrated work, before joining the present 
institute. This program was reported in [29,30]. A brief 
description of the theory used in this program is presented 
here for the sake of completeness. A modification developed 
during this work over the existing method to consider phase 
change effectively, is also explained in Section 5.4.3. 

This program has been used for computing temperature 
distribution during longitudinal welding of cylinder, as 
described in Chapter 6. 

5.2 Theory of Two Dimensional Linea r Heat Transfer 

The governing differential equation for two dimen 
sional time dependent heat conduction is given below. For 
an axisymmetric geometry the equation gets slightly modified 
and is also given below. 

For 2-D Case 

§x (K X "ai 5 + fv ( ** + ° pc B t 


(5.1) 



Boundary Conditions 
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Constant temperature T = t(x, y, t) on si, t > 


0 


(5.2) 


Convective heat loss 


9T 


Be n x + *Sr 'fl *V + 5 + h(T - T„) 


0 on s 2 , t > 0 


(5.3) 


Initial temperature 


T = T (X, Y) in D y t = 0 


(5.4) 


For Axi symmetric Case 


3T. 


3 gr (K r r "fr J + 32 (K Z r 3Z } + r Q 


= pc 


3 T 
3 t 


(5.5) 


Boundary Conditions 

Constant temperature T = T(r,z,t) on SI, t > 0 


(5.6) 


Convective heat loss 


K r n + K™ r -f£ n„ + r h (T - T ) 

r ar r Z 3Z Z 


0 on s , t > 0 


Initial temperature 


T = T (•’■*' , z) in D, t = 0 


(5.7) 


(5.8) 


The tertperature distribution over an element 'e 1 is 
assumed as 

(5.9) 


T (e) (X, Y, t) 


S Cx, Y) T . (t) = [N] PJ (ne) 


i=l 


substituting Equation (5.9) in dif 4 erential Equations 
(5.1) or (5.5) and using Galerkin' s method, we get 
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[Kj e rr } e = iQ] e - { q } e ™[K h ] e i T } e +£ K T _ } e .[K c ] { -|f } e 

(5.10) 


CD 


where 
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'x ax 3Y y ax aY- ; 


1 dX dY 


(5.11) 


Q A = / Q N ± dX dY 

D e 


(5.12) 


q, 


f q N, d S 

b 2 


(5.13) 


K h. 


/ h N. N . d X! 


ij 


S 
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i j 


(5.14) 


Kt 


go 


K 


*ij 


/ hT_N.dE 

... OD ± 

S 

S 2 


/ pCN. N . dX dY 

D e 


(5.15) 


(5.16) 


combining [Kj e and£ KjJ® and all the load vectors, 
the equation appears as 

[K] (e) ST} (ne) + [cf - 5 E '} (ne> 

The elements are assembled into global matrices and 
the final system of equations is obtained as 

_ _ „ »T, _ nji (5.17) 

[K] {T} + [C] {"ft} " 11 

, . _ f or the axisymmetric geometry 

This procedure is same tor m - 

except that 2n r dr dz replaces dX dY ana 2jl r dj r p 

d S . The final set of equations have the same form. 




The set of differential equations for time dependent 
heat conduction are given by Equation ( 5 . 17 /). 

Two procedures for solution of these equations are 
popular , namely direct integration and modal superposition. 
Since modal superposition method is not useful for nonlinear 
heat transfer problems/ direct integration methods are 

discussed here. 

In direct integration, time dependent heat conduction 
equations are integrated using a numerical step by step 
procedure without applying any transformation prior to 
integration. The heat balance is carried out only at dis- 
crete time intervals At apart over a time interval At, a 
certain variation of is assumed. 

Here finite difference approach is discussed in 
details since this has been used in the present work. 


5.3.1 Finite Difference 

The time derivative of temperature is approximated 
using finite difference. It is assumed constant over the 
time step. Let the superscripts i and i +1 denote beginning 
and end of the time step. Depending on what constant value 
is assumed, throe schemes are possible, which are listed 


as follows: 



Scheme 

Constant 

value of {||} 

1 

over 

1 . Euler' s/Explicit scheme 

t-as , 1 

1 at i 


(5.18) 

2 . Implicit scheme 



(5.19) 

3. Crank-Nicholson' s scheme 


♦ " ] 

(5.20) 


Of the three/ only Crank-Nicholson’ s scheme Is 
described here because it allows for larger time steps and 
produce;' stable oscillations in case the stability criterion 
is violated [ 13 ] . In this scheme 

ST} 14-1 = JT} 1 + -^ (5.: 

Using Equation (5.21) in Equation (5.17), we have 
after simplifications: 

[[c] + 4j~ [K]] CT} 1+1 = [ [c] - -^r [K] ]{T} X 



+ \ 

[{R} 1 + JR } 1+1 ] 

(5.22) 

or O] 

£T} ±+1 = [B] {T} 1 ■ 

+ JR} 

(5.23) 

where 




[A] 

+ 

i — i 

U 

i — i 

if 


(5.24) 

[B] 

* [C] -^r [ K ] 1 

JR} 1 = At JF} 1 

(5.25) 

and j r } 

= \ J {R} 1 + €R } i+1 

} 

(5.26) 


These equations form the time stepping algorithm. 
There is a modified version of Crank-Nicholson s 
scheme called modified Crank-Nicholson' s method [13] . If 
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average temperature of a time interval i s defined as 



!T a } = HTJ 1 4.(TJ i+1 )/2 

(5.27.) 

then 

;T} 1+1 = 2{T a } - f tj A 

(5.28) 


So above equation can be written as 


[[C] 4 -ft [K]] f2IT a } - (T)*, = [[C] - -At [K]] ( TJ i 



+ \ HR} 1 + { R} i+1 } 

(5.29) 

or 

[-ft [c] + [K]]CT a -l + {R f + [c] CT } 1 

(5,30) 


Once {T a } is calculated for every step, { T}^ + ^ is 
calculated by 


i+1 i 

{Tj = 2 {T a } - {T} (5.31) 

Present program WELTEM for nonlinear heat transfer, 
uses this solution procedure. 

5 . 4 Nonlinearity in Heat Transfer 

5.4.1 Change of Material Properties with Temperatur e 

Material properties such as heat capacitance, 
thermal conductivity, density etc. are in general temper- 
ature dependent. Many a vi.mes their temperature dependency 
can be neglected. But many mechanical processes such as 
welding, rolling etc. involve large i. nge of terrperature 
change and hence one cannot neglect change material 
properties with temperature. 
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The easiest way of considering property changes 
with temperature is, by supplying values of properties at 
different temperatures. During transient calculation one 
can linearly interpolate values of properties for a particular 
element depending upon its average temperature, before every 
time step and compute its capacitance and conductivity mat- 
rices . 

Interpolation of properties may be done depending 
upon the temperature of every gauss point or depending upon 
the average temperature of element. Former gives better 
result in terms of accuracy but involves more amount of compu- 
tations . 

5.4.2 Consideration of Thermal Radiation 

Heat loss due to thermal radiation is very signi- 
ficant when high temperature boundaries are dealt with, such 
as during welding, casting etc. This is essentially a non- 
linear problem, because of explicit temperature dependency 
of heat loss on boundary temperatures. 

For taking into consideration heat loss due to 
radiation one can define a factor h r (coefficient of radi- 
ation heat transfer) similar to h c (coefficient of convec- 
tive heat transfer) such that heat loss due to radiation is 
equal to [14] 

h r (T * T at ) = 6 r ' T at ) 

where e r is emissivity , ° r i s Stefan-Boltzmann constant 
T is the boundary temperature and T at is the atmospheric 


(5.32) 



temper <at.nro, n^rw-o 


h 


! r °r (T + T at > <T + T at> 


(5.33) 


To evaluate value of h^. , if one assumes linear 
temperature distribution within the interval t - At to 

t + At, then 

h r * 2 At ^__ At E r °r (T + T at } (T + T at } dt 


= e r a r [( T (t) 2 +T 2 t )(T(t) + T at ) + (T(t) -T(t-At)) 2 

(T(t) + iT at )] 

(5.34) 

Here T at is assumed to be constant,. Since right 
hand side of Equation (5.34) involves only known values of 
temperature, hence it can be calculated at every time step. 

is a material properties and hence may be a function of 
temperature also. 

5.4.3 Phase Change 

In different mechanical processes such as welding, 
casting etc., materials change phase during heating and 
cooling. This in turn affects the temperature distribution 

calculations during change of phase. 

The mathematical model of liquidif ication or 
solidification alloys is complicated by the fact that there 
is a moving boundary, which is the interface between th 
solid phase and the liquid phase. In pure metals th p 
change occurs at one freezing temperature, in Y 
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however, the phase change occurs 

etween two temperature 

levels called the solidus and the n 

ana the liquidus temperatures. 

Above the liquidus the alloy is oomni^i. i • 

Y 1S completely liquid, below the 

solidus it is completely solid, and in between the alloy 
consists of a mixture of solid and liquid. Hence, two 
moving boundaries, one at solidus temperature and the other 
at the liquidus temperature, exist during the solidification 
of an alloy. Keeping track of these boundaries in the finite 
element model, is very cumbersome. Therefore, in this work 
the concept of equivalent specific heat is used to approxi- 
mate the solidification process. 

This concept is from the facts that during phase 
change alloys either release or absorb large amount of heat 
without appreciable temperature change. This phenomena can 
be viewed in other way that during phase change apparent 
specific heat of alloys is very high and hence even for 
small temperature change large amount of heat is released or 
absorbed whatever the case may be. This apparent high 
specific heat C po is found as follows. 

As shown in Figure 5.1, C ps is the specific heat at 
solidus temperature and C pl is specific heat at liquidus 
temperature. is latent heat of solidification and 

and T are liquidus and solidus temperatures respectively . 

s 

Then c p can be approximated by 


C Ps + C P1 


4 


"H 


!T i - V 


(5.35) 


can 


Pe “ 2 

The equivalent specific heat value calculated above 
be used whenever the average temperature of the element 



is- between the liquidus and solidus temperature By this 
means, the effect of latent heat released during the alloy 

solidification is taken into account, and the overall heat 
balance is maintained. 

But in this method, the greatest drawback is that, 
there is a tendency of either jumping the melting temperature 
range or entering the range leaving a gap between solidus 
temperature and average temperature of the element, as shown 
in Figure 5.2, and hence not preserving full latent heat, 
associated to that element. 

Different methods are suggested to consider the phase 
change, to avoid upper mentioned drawbacks [22,40], But all 
these methods are iterative and hence requires more computer 
time. To avoid this and also to make use of upper simple 
method, during the course of the present work, a modification 
was implemented in WELTEM, which was found very effective. 

This modification is explained as follows for diff- 
erent cases: 

Different notations used are: 

T - temperature at the end of current iteration 
T^ - an estimated temperature for the next iteration, 

considering heat capacity 'pc ' corresponds to f, shows 
a change of region. It is estimated from T by adding 
temperature increment of the previous iteration. 

T£ - an improved estimation to temperature for next iter- 
ation, considering a better approximation to heat 
capacity (Pc) av * 





Case 1 


lie 


as ~pc~) 


T 1 > T s > T ( Ref * : Figure 5.3(a)) 

Heat absorbed from T f n t* /_ ^ • - . 

xrom i 1 to T. (considering heat capacity 


pc (T x - T) (5.36) 

Let this heat to the element undergoing change of 
region be a constant. This assumption is justified as very 
few elements undergo a change of region in one iteration. 

Since "2 a better estimation to , hence by 

heat balance, we have 


PC - T) = pc (T - T) + (pc) (T 0 - T ) 

o S Z S 


° r T 2 = T s + XpW <T 1 - V (5.37) 

0 

Equation (5.37) is independent of T explicitly. 

Hence defining (Pc) av as 


or 


(pc) av (T 2 - T) = pc (T s - T) + (pc) e (T 2 

Pc (T- - T) 

( PC) = — Z 

av (T 2 - T) 


T 

s 


) 


(5.38) 


Hence a better approximation to pc is ( P c ) av ,/ 


calculated by Equation (5*38) . 

For other cases, proceeding in a similar way, the 
expressions for T 2 and (p c) av ate determined as follows. 

case 2 - > Tj > T (Ref.: Figure 5.3(b)) 


and 


(P c) 


(T X - 


T l) 


( pc). 


(pc) 


av 


pc (T x - T) 

(t 2 


(5.39) 

(5.40) 



Case 3 




(5.41) 

(5.42) 

(5.43) 

(5.44) 


Any of the one pair of equations from Equations 

(5.37/) to (5.44) can be used to calculate and (pc) / 

2 av 

depending upon the case. 



PART II 


APPLICATIONS 
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CHAPTER 6 

ANALYSIS OF WELDED STRUCTURES FOR RESIDUAL STRESSES 

6 . 1 Introduction 

6 . 2 Thermal Model 

6.2.1 Quasi- St at ionary Theory 

6.2.2 Heat Input and Distribution Calculations 

6.3 Stress and Distortion Model 

6.3.1 2-D Simulation of Welding Process for Stress 
Calculations 

6.3.2 A Model for Simulating Melted Region 

6.3.3 Stress Calculations by Finite Element Method 

6.4 Solved Examples 

6.4.1 A Note on Material Properties 

6.4.2 Analysis of Butt weld of Two Plates 

6.4.3 Analysis of Longitudinal welding of 

* Cylinder 



6.1 Introduction 
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The welding of metal structures is aimed at provi- 
ding a means of joining together a number of components in 
such a way as to minimise the impairment of the properties 
of those components. Considerable effort has been expended 
on developing effective welding techniques for a large number 
of metals and alloys. Concurrent with these developments 
have been efforts aimed at identifying the various problems 
that results from welding processes, determining the strength 
of welded joints and structures subject to specified loading 
conditions, and establishing guidelines and criteria for 
most effective joint design. The information accumulated in 
these areas over the years has been overwhelmingly experi- 
mental . More recently, efforts have been directed toward 
developing analytical models to predict the thermal-mecha- 
nical response of welded structures* 

Analytical formulations for predicting temperature 
distribution during welding are frequently obtained by ass- 
uming that the thermal energy supplied by the weld process 
can be idealized as a point or line heat source. Early work 
in this area focussed on quasi-stationary , transient temper- 
ature distributions resulting from a point heat source trave- 
lling at a constant speed along a line on an infinitely 
plate. The line heat source and moving point heat source 
have been used to obtain good correlation with temperature 
measurements for welding conditions. As an alternative to 
the analytical approach, numerical techniques like finite 
difference and finite element methods have also been used. 



Analytical models for calculating stresses during 

welding have been developed based on the line heat source 

and moving point heat source Th^cA +.^ . . 

source. These techniques generally 

show good agreement between residual 0 +- v 

residual stresses of the model 

and experimental data but are often limited to single pass 
welds. The complexities inherent in performing analytical 


studios of welding stresses and distortions for various weld 
configurations - in particular, irregular geometries, in- 
elastic material response and material loading and reloading 
- suggest the use of numerical methods to-calculate, residual 
stresses and distortions. Because of its generality, finite 
element method has been applied to study the behavior of the 


welded structures. 

To show analysis capability .of program THESIS, which 
has been developed during the present work, depending upon 
the theory presented In Chapters 3 and 4, stress transients 
during butt welding of two plates and longitudinal welding 
of a cylinder are presented here. In fact programs WELTEM 
and THESIS combined provide a complete set to compute temp- 
erature and stress transients during welding, casting, 
annealing etc . 


6.2 Thermal Model 

6.2.1 Qua si- Stationary Theory 

Calculation of the transient temperature distribu- 
tion is based on the attainment of quasi-stationary condi- 
tions , which are developed when the welding heat source is 
moving at constant speed on a regular path, and also end 




effects resulting from either initiation or ter mination of 
the heat source are neglected. The temperature distribution 
is then stationary with respect to a moving coordinate system 

whose origin coincides with the point of application of the 
heat source . 

Considering the planar weld shown in Figure fe.l. The 
temperature at any point in the weldment is expressed funct- 
ionally as : 

T(x, y , z, t) = T(x, y, z, -vt) (6.1) 

where v is the welding speed. Thus, given the transient 
temperature distribution at any one section of the weldment 
defined, say, by z = 0, the tenperature at any other section 
is determined by an appropriate shift of the time scale as 
follows : 

T (x, y, z, t) = T(x, y, o, t - z/v.) (6.2) 

The problem is therefore,, reduced to finding the 
two dimensional, unsteady temperature field at a section 
normal to the weld line. A planar analysis may be used for 
this purpose when the weld speed, relative to a character- 
istics diffusion rate for the material, is sufficiently high 
so that the amount of heat conducted ahead of the weld torch 
is very small relative to the total heat input. In this 
case, the net heat flow across any infinitesimally thin 
slice of the weldment normal to the weld line is assumed to 

be negligible relative to the heat being diffused within 

3 3T 

the slice itself, that is, the term ^ ±S 

neglected in the heat conduction equation* Two dimensional 
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thermal analysis at the section z = 0, normal to the direc- 
tion of welding, is thus treated in the present work. 

6.2.2 Hggt Input and Distribution Cal cul ations 

Perhaps the most critical input data required for 
welding thermal analysis are the parameters necessary to 
describe the heat input to the weldment from the arc. Not 
only the magnitude, but the distribution of heat input will 
influence the dimensions of the weld metal and heat-affected 
zones, the cooling rates, and the peak temperature distri- 
bution, as well as the temperature gradients necessary to 
calculate stresses and distortions. The welding parameters 
required to formulate the heat flux boundary conditions are 
the magnitude of heat input from the arc Q, the- distribution 
of the heat input, characterized by a length parameter , 
and the weld speed v. Introducing the arc ef f iciency V , Q 
is found simply from the formula 

q = r? s I (6.3) 

where E and I are the arc voltage and current, respectively. 
The heat deposited per unit length of weld is merely Q/v. 

It must be recognized that a priori determination of the 
magnitude and distribution of heat input cannot, in general, 
be made due to the lack of knowledge regarding energy trans- 
fer from the arc to the workpiece. Investigations of the 
physics of the welding arc are required to shed light on 

this area. 

The heat from the welding arc is, at any given time, 

, . , , . , _ j mirfacB of the weldment as 3. 

assumed to be deposited on tne surrey 



radially symmetric normal distribution function [ll] . 

Letting r be the distance from the center of the heat source 
which is coincident with the axis of the electrode, the heat 
flux q is given by 


= % e '° ^ ( 6 . 4 ) 

where q^ and c are constants determined by the magnitude 
and distribution of the heat input. 

The heat input parameters Q and r are defined by [ll] 

OD 

Q = 2 7t f q(r) r dr ( 6 . 5 ) 

o 


and q(r) = 0.05 < 1 0 

Then Equation (6.4) becomes 


( 6 . 6 ) 


m ~[>( £ ) ] 

q(r) = — 3 e ^ (6.7) 

nr 

Here r defines the region in which 95 percent of the 
heat flux is deposited. Let x be the distance from the weld 
line in the section Z = 0 and t = 0 is the time at which the 


centre of the heat source (electrode) passes over this 
section. The heat flux distribution on the surface of the 
weldments, directed parallel to the electrode, is given by 

[u] 


-[3 0 ] -[ 3(0 ] 


q(x, t) 


M- e f 


nr 


( 6 . 8 ) 


The three welding parameters Q, r and v are all 


embodied in this formula. 



6-3 .Stress and Dist ort j on Model 


6 - 3 - 1 

A complete thorough determination of the mechanical 
response due to the welding thermal cycle, as determined by 
the quasi-stationary temperature analysis, would require a 
full three dimensional incremental plasticity analysis with, 
at the least, recalculation of the stiffness coefficients at 
each time step . Computer run times would be very large and 
the cost, therefore, will be prohibitive. Hence in the 
present work, the section normal to the weld direction, is 
analysed. This implies that all sections normal to the weld 
line remain plane during the entire welding process. Though 
this assumption is probably adequate in regions somewhat 
removed from the weld puddle, it may not be valid in the 
neighborhood of the molten metal. The approach employed at 
present therefore prevents calculation of distortions ahead 
of tha welding arc,, and may inhibit accurate computations of 
deformations in the immediate vicinity of the weld puddle. 

It nevertheless enables the essential features of the 
mechanical response during cool down to be modelled, and the 
resultant residual stresses and distortions to be calculated. 

6*3.2 A Model for Simulating Melted Domain 

The stiffness becomes negligible as the melting 
range is approached, but the material also is relatively 
incompressible in the molten state. In order to model this 
behavior and, at the same time, to avoid possible ill-behavior 



ol the solution process and thus facilitate convergence, 
the elastic modulus approaches zero and the poisson' s ratio 
one half , in such a way that the bulk modulus is maintained 
at its room temperature value through the entire high temper- 

ature range, including the phase transition and the liquid 
state. 

For illustration, since for plane str ai n case stress- 
strain matrix [d] is given by 


[D] 


E 


(1 + V ) (1 - 2V) 


l-v 

V 

o 


V 


0 

1-2V 


for melted region enforcing the conditions E 

E 


0, V -0.5 


such that *r 


i d - m * v Hence 


[f] 




i 

o 


i 

i 

0 


0 

0 

0 


(6.9) 


Since compressibility has no effect on the general- 
ized plastic strain, no plastic strain is, as expected, 
accumulated in the weld puddle. 


6.3.3 Stress Calculations b y Finite Element Metho d 

Stress transient computations during welding of 
a plane stress or plane strain or axi symmetric structure 
can be done by conputer program THESIS . Theory of isopara- 
metric finite element is given in Appendix 1 and calcul 
ations of material matrix and load vector are shown for a 
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2-D thermoplastic structure in Chapter 3. Solution proce- 
dure along with elastic-plastic transient calculations 
adopted by program THESIS are also shown in Chapter 3. 

6 „ 4 S olved Examples 

6.4* 1 A Note on Material Properties 

The weldment material properties employed are those 
as used by Friedman [20], which are essentially for Inconel 
alloy 600, an alloy of nickel, chromium and iron. 

V ariation of Material Properties with Temperature - Varia- 

Th 

tions of E,V t O , c , K, e with temperature are shown in 

mk Jr 

Figure 6.2, taken from [20]. 

3 

De nsity and Latent Heat - Density is taken as 8430 Kg/m 
and latent heat as 309 kJ/Kg. The solidus and liquidus 
temperatures are taken 1630°K and 1690°K respectively. 


Harde n ing Rule - Isotropic strain hardening of the yield 
surface is assumed, with the rule 


1 _ 

"c? = o y (s PL /0.002) 11 


( 6 . 10 ) 


Effect of Recrystallization - Following empirical expre- 
ssion is used to approximate the reduction of previously 
accumulated plastic strain at high temperatures 


PL 

ij 


f(T) e 


PL 

ij 


( 6 . 11 ) 


where 


f (T) 


0.02 (T-T m ) 


if T < T 


m 


0 


if T > T, 


m 


( 6 . 12 ) 
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where = melting temperature = 


131 

1630°K. 

6 * ' 4 ' 2 •^^s^o£_.But ^Weld of two Plates 

Aim to analyse this case here is to compare the 

results of THESIS and by other authors [20] using different 
analysis procedure. 

Problem consists of a 2.54 mm thick flat plate of 
width 150 mm. The section is assumed in a plane strain 
condition, with the heat input data as (Ref.: Equation (6.8)) 

Q = 703 W 

r = 5.08 mm 

v = 2.12 mm/sec 

Other material properties are taken, as explained in 
Section 6.4.1. Temperature transient analysis of this case 
was done by WELTEM and was reported in [30] by the present 
author. Since object of the analysis of this case, during 
the present work, is to test the formulation given in 
Chapter 3, hence temperature distribution was used from 
Figure 4 of [20] , instead of solving by WELTEM, to avoid any 
incompatibility in temperature distribution during the 
analysis . 

In the present analysis only 50 mm from the line 
of welding is considered for finite element discretization, 
as no stresses are developed after this length, as shown in 
Figure 7. of [20]. There are 30 elements and 125 nodes with 
two degrees of freedom per node used for the finite element 
model. Figure 6.3 shows this discretization as plotted by 


computer. 
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Figures 6.4 to 6.6 show the mid surface z direction 
stress variation at times -0.016 min, 0.15 min and 0.3 min 
respectively. Comparison with [20] shows, a close agreement, 
though the formulations used in both cases are entirely diff- 
erent. A careful comparison reveals that stresses calculated 
by THESIS in heat affected zone are lower than shown in [20], 
This difference of stresses is reported, without giving any 
proof, by Rybick.i ot al in their paper 'A Finite Element 
Model for Residual Stresses in Girth Butt Welded Pipes' 
published in [7]. Figure 6.7 shows complete mid surface o 
variations at different timings . 

6.4.3 Analysis of Longitudinal Welding of a Cylinder 

Program WELTEM and THESIS have been used to analyse 
longitudinal welding of a cylinder for temperature and stress 

transients . 

For this, a cylinder of 50 mm radius is discretized 
with 50 finite elements and 205 nodes. Material properties 
and heat input data are again taken as given in Section 6.4.1. 

Figure 6.8 shows the finite element model of the 
cross-section of cylinder, assumed in a plane strain condi- 
tion. Figures 6.9 and 6.10 show temperature transient and 
mid surface stress ^ variation along the circumference for 

different time intervals. 

Stress transients show that upon initial heat up, 
the localization of severe temperature gradients in the 
immediate vicinity of the weld line produces compressive 
Yielding in this region. A temperature increases, the 
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yiold strength quickly decreases untill, at the melting 
point, it in negligible. During the time period prior to 
solidification of the weld metal, all material outside the 
puddle is in compression with that region immediately adja- 
cent to the molten zone in a plastic stats of stress. Upon 
sol ic Jit leal, ion, tho fusion zons material yislds initially 
in compression at rather low stress lsvsls « Upon ftrther 
cooling tonal le stresses are produced after unloading at 
heat affected zone. Hence, as cooldown proceeds, that 
portion of the weldment that is In tension longitudinally 
grows steadily until, at final cool-down, the residual 
longitudinal, stresses, which are appreciable only in the 
region within about 80 mm of the weld centerline, are 
completely tensile. This is due to the plane strain rest- 
rictions placed on the analysis. The residual stresses in 
the plastic region exceed the room temperature yield stren- 
gth because of material strain hardening, which is such 
that the initial short time compressive yielding produces 
an expansion of the yield surface followed, during cool- 
down, by yielding in tension at a stress level higher than 
the room temperature yield strength. 



CHAPT ER 7 
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J! 



^g^ENg-OF WELDED STRTTrvmTp^ 


7 .1 introduction 


7.2 Annealing of: Butt Welded Plate 

t 

7.3 Annealing of Butt Welded Cylinder 





7.1 Introduction 


■1.L has bean shown in Chapter 6, that the residual 
stresses at the welded joint are generally above the yield 
point at room temperature. These high residual stresses 
reduce, the life of the component. Therefore it is necessary 
to rc drive these stresses either by heat treatment or by 
some othtir mechanical method such as vibratory stress relief. 
Annea l i.ng Ah a type of popular heat treatment process y in 
which welded structure is heated below recrystallisation 
temperature slowly and is kept for a long time at that 
temperature and then cooled down to room temperature slowly. 

Relaxation of residual stresses occurs initially 
during heating due to the decrease of yield point and Young* s 
modulus and then due to generation of creep strain. 

The program THESIS can be used to predict stress 
transient in the welded structure during annealing . To show 
the capability of THESIS , results of annealing of two struc- 
tures arc presented here. In both cases creep properties 
are taken as described in Section 4.5.1. 




Butt Welded Plate 


Butt weld of two plates, as described in Section 
6. 4. ,2, is solved for annealing. The plates were assumed to 
be heated up to 600° C at the rate of 600° C/hr and held at 
this temperature for 2 hrs and then cooled to room temper- 
ature at the rate of 600° C/hr. 

The residual stress obtained from the analysis of 
Section 6.4.2 was used as input for annealing analysis. 






STRESS m Mh/ ytv 
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stress transients during complete annealing at different 
time interval tt are shown in Figure 7.1, as plotted by computer 

7 . 3 Annealing of Butt Welded Cylinder 

The example of longitudinal weld of a cylinder, 
described in Section 6.4.3, was then analysed for annealing. 
For thin, again the cylinder is heated to a 600° c temperature 
in one ho vu and then held at this temperature for two hours 
and is th* *i i cooled to room temperature in one hour time. The 
stress transients during this annealing cycle is plotted by 
computer for different time intervals are shown in Figure 7.2. 

In both the above analysis, it is seen that residual 
stress decreases slowly till 400° C and decreases rapidly 
afterwards because of the generation of creep strain. No 
significant increase in residual stress is obtained while 
cooling the structure, though yield point and Young' s modulus 
increase. This is because of the generation of large creep 
strain during the two hours of holding time at constant 
temperature which completely dominates the process of stress 
relief . 



CHAPTER 8 


CONCLUSIONS AND FIJTTTPg 


Conclusions 


Future Work 



B.l Conclus ions 

A Ww«»-pU*tio formulation with creep for plane 
strain case is derived and used to solve different problems. 

along with welding and annealing. 

Finite element technique is used to solve the whole 
stru duro *nd found to be very useful for nonlinear problems 
also T,u «3*«nb stiffnima solution procedure is used for time 
stepping algorithm. It is seen that drift from yield surface 
can be mad© very small provided load increment in every 
iteration is kept low and remaining unbalance force is added 
to the load vector of the next iteration. This way internal 
iterations can be avoided and hence less computer time. 

Since in plane-strain formulation, stiffness as well as load 
vector are functions of stress field of the previous itera- 
tion. it is necessary to compute stress field very accurately 
as small error diverges through further iterations. This 
precaution is more valid whenever there is elastic-plastic 
or plastic- elastic transition. Hence internal iterations 
should be done to bring the stress of a transition point to 
the yield surface. 

Though the thermo-plastic formulation and solution 
procedure for the butt welding of two plates adopted m 
program THESIS and in [20] are entirely different, a close 
agreement of the temporal stresses distribution computed by 

both methods is obtained. Stresses at the 

zone calculated of 'THESIS is found lower than reported in 
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StrGSS transients during annealing of butt welded 
plates oi cylinder show slow decrease of residual stresses 
during initial heat up till 400° c. After this stress© ~ 
decreaso rapidly duo to the generation of creep strain. 

When the structure is kept at 600° c for two hours and then 
cooled to room temperature, a negligible increase in stresses 
occur . this is due to the generation of large creep strain 
during holding time, which dominates the while stress 
distribution* 

8.2 Future Work 

Present computer program THESIS can be modified for 
taking into consideration Multipass welding of structures. 

Also since during the entire range of temperature , there is 
phase transformation of iron, coefficient of thermal expan- 
sion changes significantly within phase transformation 
temperature range. Present computer program can also be 
modified to take into account change in coefficient of thermal 

expansion due to this factor. 

Experimental verifications of residual stress distri- 
bution during butt welding should be done to test validity 
of the theory and solution procedure used. Though some 
efforts have been put by different authors [10,11,19,21] to 
measure residual stresses during butt welding, but present 
author intends to do experiments for electron beam welding. 
Because of less radiation heat loss, good control on he 
input data, less heat affected zone, a better mathematical 
simulation is obtained for electron beam welding 



a better comparison of stress and temperature variations is 

expect* *d . 

Thou* jh present computer program has only been used 
to solve welding and annealing problems, stress transients 
during casting can also be computed by this program, with 
little modif ications/ capability of the present program can 
be extended to predict stress transients during other mecha- 
nical processus such as grinding etc. 
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APPENDIX 1 

ISOPARAMETRIC FINITE ELEMENT 2-D FORMULATIONS 


Element stiffness matrix of a finite element is given 

by 

[K] <e) = l of [D] [b] dv 

V 

where [b] is called strain-displacement matrix and [d] is 
called stress-strain matrix. 

Strains at a point in finite element is given by 

{&} = [b] {d} 


where { e } for plane stress case and plane strain case, is 

given by 


. W = U x s y 

and for axisymmetric case 



UJ = Le x e y e xy e 0 ^ 

{d} is the vector of displacements at a point, hence {d} = < 
So matrix [5] for plane stress case and plane strain case is 
given by 
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and for axisymmetric case 


where r is the radius of the point under consideration. 
In case of isoparametric formulation 


S N, u, ; V 
i=l 


S N . v. 
1=1 


N i X i 


and y ■ s N i y i 
i=l 

where n is the number of nodes per element and N ± is the no 
shape function . Also v v., x,, y, are the corresponding 

nodal values. 

Using these interpolations for equations for 
strains, we have for plane stress and plane strain case 


[B] {53 
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Also for axi symmetric caso 
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Also since 
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where £ and V are the local coordinates used to define 
shape functions , we have 


f . -» 

3 N, 


r"* 

3N. 

| 

a 1 


ax ay 

t 

1 

1 

a i 


si 3 e 


ax 

1 

% /V KT 

> 


-<! 

KT 

v = t j ] < 

d am 

i 

I 

$x ay 


9 i 


3 H 


a r? ’ 3 ti 


ay 



iNfj 

X 


ax 
3N . 


i 


ay 


where [j] is called Jacobian matrix. Hence 


a n ~ 
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f8 N l1 

3x | 


3 1 
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3 N , | 
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an . 


Using interpolation defined earlier for coordinate 


at a point in terms of nodal coordinates, we have 
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calculated by upper expressions shown. 

Matrix [d] for different cases are given by for plane 


stress case; 


1 - v 


for plane strain case: 


+ V ) (1 - 2V 


V ' l-y 0 

0 0 ~r 


1 - 2 * 

2 


for axi symmetric case: 
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